diff --git a/doc/src/Fortran.rst b/doc/src/Fortran.rst index ec54a594d1..4c450a0c08 100644 --- a/doc/src/Fortran.rst +++ b/doc/src/Fortran.rst @@ -30,7 +30,7 @@ as a static or dynamic 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``, ``mpiexec`` or equivalent. Please also note that the order +``mpirun``, ``mpiexec``, 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 :f:mod:`LIBLAMMPS` module that is imported by the Fortran code that uses the interface. A working example @@ -39,21 +39,12 @@ can be found together with equivalent examples in C and C++ in the .. versionadded:: 9Oct2020 -.. admonition:: Work in Progress - :class: note - - This Fortran module is work in progress and only the documented - functionality is currently available. The final implementation should - cover the entire range of functionality available in the C and - Python library interfaces. - .. note:: - A contributed (and more complete!) Fortran interface that more - closely resembles the C library interface is available in the - ``examples/COUPLE/fortran2`` folder. Please see the ``README`` file - in that folder for more information about it and how to contact its - author and maintainer. + A contributed Fortran interface that more closely resembles the C library + interface is available in the ``examples/COUPLE/fortran2`` folder. Please + see the ``README`` file in that folder for more information about it and how + to contact its author and maintainer. ---------- @@ -79,7 +70,7 @@ function and triggered with the optional logical argument set to USE LIBLAMMPS ! include the LAMMPS library interface IMPLICIT NONE TYPE(lammps) :: lmp ! derived type to hold LAMMPS instance - CHARACTER(LEN=*), PARAMETER :: args(3) = & + CHARACTER(LEN=12), PARAMETER :: args(3) = & [ CHARACTER(LEN=12) :: 'liblammps', '-log', 'none' ] ! create a LAMMPS instance (and initialize MPI) @@ -136,11 +127,11 @@ 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 +:f:func:`commands_string`) have to make copies of the strings passed as arguments so that they can be modified to be compatible with the requirements of strings in C without affecting the original strings. Those copies are automatically deleted after the functions return. -Below is a small demonstration of the uses of the different functions: +Below is a small demonstration of the uses of the different functions. .. code-block:: fortran @@ -286,6 +277,8 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. :ftype extract_fix: function :f extract_variable: :f:func:`extract_variable` :ftype extract_variable: function + :f set_variable: :f:subr:`set_variable` + :ftype set_variable: subroutine :f gather_atoms: :f:subr:`gather_atoms` :ftype gather_atoms: subroutine :f gather_atoms_concat: :f:subr:`gather_atoms_concat` @@ -298,6 +291,16 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. :ftype scatter_atoms_subset: subroutine :f gather_bonds: :f:subr:`gather_bonds` :ftype gather_bonds: subroutine + :f gather: :f:subr:`gather` + :ftype gather: subroutine + :f gather_concat: :f:subr:`gather_concat` + :ftype gather_concat: subroutine + :f gather_subset: :f:subr:`gather_subset` + :ftype gather_subset: subroutine + :f scatter: :f:subr:`scatter` + :ftype scatter: subroutine + :f scatter_subset: :f:subr:`scatter_subset` + :ftype scatter_subset: subroutine :f create_atoms: :f:subr:`create_atoms` :ftype create_atoms: subroutine :f find_pair_neighlist: :f:func:`find_pair_neighlist` @@ -308,6 +311,8 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. :ftype find_compute_neighlist: function :f neighlist_num_elements: :f:func:`neighlist_num_elements` :ftype neighlist_num_elements: function + :f neighlist_element_neighbors: :f:subr:`neighlist_element_neighbors` + :ftype neighlist_element_neighbors: subroutine :f version: :f:func:`version` :ftype version: function :f get_os_info: :f:subr:`get_os_info` @@ -358,6 +363,22 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. :ftype encode_image_flags: function :f decode_image_flags: :f:subr:`decode_image_flags` :ftype decode_image_flags: subroutine + :f set_fix_external_callback: :f:subr:`set_fix_external_callback` + :ftype set_fix_external_callback: subroutine + :f fix_external_get_force: :f:func:`fix_external_get_force` + :ftype fix_external_get_force: function + :f fix_external_set_energy_global: :f:subr:`fix_external_set_energy_global` + :ftype fix_external_set_energy_global: subroutine + :f fix_external_set_virial_global: :f:subr:`fix_external_set_virial_global` + :ftype fix_external_set_virial_global: subroutine + :f fix_external_set_energy_peratom: :f:subr:`fix_external_set_energy_peratom` + :ftype fix_external_set_energy_peratom: subroutine + :f fix_external_set_virial_peratom: :f:subr:`fix_external_set_virial_peratom` + :ftype fix_external_set_virial_peratom: subroutine + :f fix_external_set_vector_length: :f:subr:`fix_external_set_vector_length` + :ftype fix_external_set_vector_length: subroutine + :f fix_external_set_vector: :f:subr:`fix_external_set_vector` + :ftype fix_external_set_vector: subroutine :f flush_buffers: :f:subr:`flush_buffers` :ftype flush_buffers: subroutine :f is_running: :f:func:`is_running` @@ -407,7 +428,7 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. USE LIBLAMMPS USE MPI_F08 TYPE(lammps) :: lmp - lmp = lammps(MPI_COMM_SELF%MPI_VAL) + lmp = lammps(comm=MPI_COMM_SELF%MPI_VAL) END PROGRAM testmpi .. f:type:: lammps_style @@ -575,9 +596,9 @@ Procedures Bound to the :f:type:`lammps` Derived Type .. code-block:: fortran ! code to start up - logical :: periodic(3) + LOGICAL :: periodic(3) ! code to initialize LAMMPS / run things / etc. - call lmp%extract_box(pflags = periodic) + CALL lmp%extract_box(pflags = periodic) -------- @@ -643,8 +664,8 @@ Procedures Bound to the :f:type:`lammps` Derived Type USE MPI_F08 USE LIBLAMMPS - TYPE (lammps) :: lmp - TYPE (MPI_Comm) :: comm + TYPE(lammps) :: lmp + TYPE(MPI_Comm) :: comm ! ... [commands to set up LAMMPS/etc.] comm%MPI_VAL = lmp%get_mpi_comm() @@ -684,7 +705,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type 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 + 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. @@ -733,7 +754,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type CHARACTER(LEN=:), ALLOCATABLE :: str lmp = lammps() CALL lmp%command('units metal') - ALLOCATE ( CHARACTER(LEN=80) :: str ) + ALLOCATE(CHARACTER(LEN=80) :: str) str = lmp%extract_global('units') str = TRIM(str) ! re-allocates to length len_trim(str) here PRINT*, LEN(str), LEN_TRIM(str) @@ -771,8 +792,8 @@ Procedures Bound to the :f:type:`lammps` Derived Type Note that this function actually does not return a pointer, but rather associates the pointer on the left side of the assignment to point to internal LAMMPS data. Pointers must be of the correct type, kind, and - rank (e.g., ``INTEGER(c_int), DIMENSION(:)`` for "type", "mask", or "tag"; - ``INTEGER(c_int64_t), DIMENSION(:)`` for "tag" if LAMMPS was compiled + rank (e.g., ``INTEGER(c_int), DIMENSION(:)`` for "type", "mask", or "id"; + ``INTEGER(c_int64_t), DIMENSION(:)`` for "id" if LAMMPS was compiled with the ``-DLAMMPS_BIGBIG`` flag; ``REAL(c_double), DIMENSION(:,:)`` for "x", "v", or "f"; and so forth). The pointer being associated with LAMMPS data is type-, kind-, and rank-checked at run-time. @@ -897,7 +918,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type TYPE(lammps) :: lmp REAL(c_double), DIMENSION(:), POINTER :: COM ! code to setup, create atoms, etc. - CALL lmp%compute('compute COM all com') + CALL lmp%command('compute COM all com') COM = lmp%extract_compute('COM', lmp%style%global, lmp%style%type) will bind the variable *COM* to the center of mass of the atoms created in @@ -1071,10 +1092,10 @@ Procedures Bound to the :f:type:`lammps` Derived Type - Per-atom array In the case of global data, this function returns a value of type - ``real(c_double)``. For per-atom or local data, this function does not + ``REAL(c_double)``. For per-atom or local data, this function does not return a value but instead associates the pointer on the left side of the assignment to point to internal LAMMPS data. Pointers must be of the correct - data type to point to said data (i.e., ``REAL(c_double)``) and have + type and kind to point to said data (i.e., ``REAL(c_double)``) and have compatible rank. The pointer being associated with LAMMPS data is type-, kind-, and rank-checked at run-time via an overloaded assignment operator. @@ -1178,8 +1199,8 @@ Procedures Bound to the :f:type:`lammps` Derived Type .. f:function:: extract_variable(name[,group]) - This function calls :cpp:func:`lammps_extract_variable` and returns a scalar, - vector, or string containing the value of the variable identified by + This function calls :cpp:func:`lammps_extract_variable` and returns a + scalar, vector, or string containing the value of the variable identified by *name*. When the variable is an *equal*-style variable (or one compatible with that style such as *internal*), the variable is evaluated and the corresponding value returned. When the variable is an *atom*-style variable, @@ -1240,6 +1261,22 @@ Procedures Bound to the :f:type:`lammps` Derived Type -------- +.. f:subroutine:: set_variable(name, str) + + Set the value of a string-style variable. + + .. versionadded:: 3Nov2022 + + This function assigns a new value from the string *str* to the string-style + variable *name*\ . If *name* does not exist or is not a string-style + variable, an error is generated. + + :p character(len=*) name: name of the variable + :p character(len=*) str: new value to assign to the variable + :to: :cpp:func:`lammps_set_variable` + +-------- + .. f:subroutine:: gather_atoms(name, count, data) This function calls :cpp:func:`lammps_gather_atoms` to gather the named @@ -1254,21 +1291,26 @@ Procedures Bound to the :f:type:`lammps` Derived Type of atoms, see :f:func:`gather_atoms_subset`. The *data* array will be ordered in groups of *count* values, sorted by atom - ID (e.g., if *name* is *x* and *count* = 3, then *data* = x[1][1], x[2][1], - x[3][1], x[1][2], x[2][2], x[3][2], x[1][3], :math:`\dots`); - *data* must be ``ALLOCATABLE`` and will be allocated to length - (*count* :math:`\times` *natoms*), as queried by - :f:func:`extract_setting`. + ID (e.g., if *name* is *x* and *count* = 3, then *data* = [*x*\ (1,1), + *x*\ (2,1), *x*\ (3,1), *x*\ (1,2), *x*\ (2,2), *x*\ (3,2), *x*\ (1,3), + :math:`\dots`]); *data* must be ``ALLOCATABLE`` and will be allocated to + length (*count* :math:`\times` *natoms*), as queried by + :f:func:`get_natoms`. + + This function is not compatible with ``-DLAMMPS_BIGBIG``. :p character(len=\*) name: desired quantity (e.g., *x* or *mask*) :p integer(c_int) count: number of per-atom values you expect per atom (e.g., 1 for *type*, *mask*, or *charge*; 3 for *x*, *v*, or *f*). Use *count* = 3 with *image* if you want a single image flag unpacked into *x*/*y*/*z* components. - :p real(c_double) data [dimension(:),allocatable]: array into which to store + :p data: array into which to store the data. Array *must* have the ``ALLOCATABLE`` attribute and be of rank 1 (i.e., ``DIMENSION(:)``). If this array is already allocated, it will be - reallocated to fit the length of the incoming data. + reallocated to fit the length of the incoming data. It should have type + ``INTEGER(c_int)`` if expecting integer data and ``REAL(c_double)`` if + expecting floating-point data. + :ptype data: polymorphic,dimension(:),allocatable :to: :cpp:func:`lammps_gather_atoms` .. note:: @@ -1308,15 +1350,20 @@ Procedures Bound to the :f:type:`lammps` Derived Type :f:func:`gather_atoms`; for a similar array but for a subset of atoms, see :f:func:`gather_atoms_subset`. + This function is not compatible with ``-DLAMMPS_BIGBIG``. + :p character(len=\*) name: desired quantity (e.g., *x* or *mask*) :p integer(c_int) count: number of per-atom values you expect per atom (e.g., 1 for *type*, *mask*, or *charge*; 3 for *x*, *v*, or *f*). Use *count* = 3 with *image* if you want a single image flag unpacked into *x*/*y*/*z* components. - :p real(c_double) data [dimension(:),allocatable]: array into which to store + :p data: array into which to store the data. Array *must* have the ``ALLOCATABLE`` attribute and be of rank 1 (i.e., ``DIMENSION(:)``). If this array is already allocated, it will be - reallocated to fit the length of the incoming data. + reallocated to fit the length of the incoming data. It should have type + ``INTEGER(c_int)`` if expecting integer data and ``REAL(c_double)`` if + expecting floating-point data. + :ptype data: polymorphic,dimension(:),allocatable :to: :cpp:func:`lammps_gather_atoms_concat` -------- @@ -1330,7 +1377,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type .. versionadded:: 3Nov2022 This subroutine gathers data for the requested atom IDs and stores them in a - one-dimensional array allocated by the user. The data will be ordered by + one-dimensional allocatable array. The data will be ordered by atom ID, but there is no requirement that the IDs be consecutive. If you wish to return a similar array for *all* the atoms, use :f:func:`gather_atoms` or :f:func:`gather_atoms_concat`. @@ -1338,9 +1385,12 @@ Procedures Bound to the :f:type:`lammps` Derived Type The *data* array will be in groups of *count* values, sorted by atom ID in the same order as the array *ids* (e.g., if *name* is *x*, *count* = 3, and *ids* is [100, 57, 210], then *data* might look like - [x(1,100), x(2,100), x(3,100), x(1,57), x(2,57), x(3,57), x(1,210), - :math:`\dots`]; *ids* must be provided by the user, and *data* must be - of rank 1 (i.e., ``DIMENSION(:)``) and have the ``ALLOCATABLE`` attribute. + [*x*\ (1,100), *x*\ (2,100), *x*\ (3,100), *x*\ (1,57), *x*\ (2,57), + *x*\ (3,57), *x*\ (1,210), :math:`\dots`]; *ids* must be provided by the + user, and *data* must be of rank 1 (i.e., ``DIMENSION(:)``) and have the + ``ALLOCATABLE`` attribute. + + This function is not compatible with ``-DLAMMPS_BIGBIG``. :p character(len=\*) name: desired quantity (e.g., *x* or *mask*) :p integer(c_int) count: number of per-atom values you expect per atom @@ -1349,10 +1399,13 @@ Procedures Bound to the :f:type:`lammps` Derived Type *x*/*y*/*z* components. :p integer(c_int) ids [dimension(:)]: atom IDs corresponding to the atoms to be gathered - :p real(c_double) data [dimension(:),allocatable]: array into which to store + :p data: array into which to store the data. Array *must* have the ``ALLOCATABLE`` attribute and be of rank 1 (i.e., ``DIMENSION(:)``). If this array is already allocated, it will be - reallocated to fit the length of the incoming data. + reallocated to fit the length of the incoming data. It should have type + ``INTEGER(c_int)`` if expecting integer data and ``REAL(c_double)`` if + expecting floating-point data. + :ptype data: polymorphic,dimension(:),allocatable :to: :cpp:func:`lammps_gather_atoms_subset` -------- @@ -1372,13 +1425,14 @@ Procedures Bound to the :f:type:`lammps` Derived Type The *data* array needs to be ordered in groups of *count* values, sorted by atom ID (e.g., if *name* is *x* and *count* = 3, then - *data* = [x(1,1) x(2,1) x(3,1) x(1,2) x(2,2) x(3,2) x(1,3) :math:`\dots`]; - *data* must be of length (*count* :math:`\times` *natoms*). + *data* = [*x*\ (1,1), *x*\ (2,1), *x*\ (3,1), *x*\ (1,2), *x*\ (2,2), + *x*\ (3,2), *x*\ (1,3), :math:`\dots`]); *data* must be of length *natoms* + or 3\*\ *natoms*. :p character(len=\*) name: quantity to be scattered (e.g., *x* or *charge*) :p data: per-atom values packed in a one-dimensional array containing the data to be scattered. This array must have length *natoms* - (e.g., for *type* or *charge*) or length *natoms*\ :math:`\times 3` + (e.g., for *type* or *charge*) or length *natoms*\ :math:`{}\times 3` (e.g., for *x* or *f*). The array *data* must be rank 1 (i.e., ``DIMENSION(:)``) and be of type ``INTEGER(c_int)`` (e.g., for *mask* or *type*) or of type ``REAL(c_double)`` (e.g., for *x* or *charge* or *f*). @@ -1403,8 +1457,9 @@ Procedures Bound to the :f:type:`lammps` Derived Type The *data* array needs to be organized in groups of 1 or 3 values, depending on which quantity is being scattered, with the groups in the same order as the array *ids*. For example, if you want *data* to be the array - [x(1,1) x(2,1) x(3,1) x(1,100) x(2,100) x(3,100) x(1,57) x(2,57) x(3,57)], - then *ids* would be [1 100 57] and *name* would be *x*. + [*x*\ (1,1), *x*\ (2,1), *x*\ (3,1), *x*\ (1,100), *x*\ (2,100), + *x*\ (3,100), *x*\ (1,57), *x*\ (2,57), *x*\ (3,57)], then *ids* would be + [1, 100, 57] and *name* would be *x*. :p character(len=\*) name: quantity to be scattered (e.g., *x* or *charge*) :p integer(c_int) ids [dimension(:)]: atom IDs corresponding to the atoms @@ -1426,7 +1481,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type .. versionadded:: 3Nov2022 - This function copies the list of all bonds into an allocated array. + This function copies the list of all bonds into an allocatable array. The array will be filled with (bond type, bond atom 1, bond atom 2) for each bond. The array is allocated to the right length (i.e., three times the number of bonds). The array *data* must be of the same type as the LAMMPS @@ -1450,10 +1505,10 @@ Procedures Bound to the :f:type:`lammps` Derived Type INTEGER :: i ! other commands to initialize LAMMPS, create bonds, etc. CALL lmp%gather_bonds(bonds) - bonds_array(1:3,1:SIZE(bonds)/3) => bonds + bonds_array(1:3, 1:SIZE(bonds)/3) => bonds DO i = 1, SIZE(bonds)/3 - WRITE(OUTPUT_UNIT,'(A,1X,I4,A,I4,1X,I4)') 'bond', bonds(1,i), & - '; type = ', bonds(2,i), bonds(3,i) + WRITE(OUTPUT_UNIT,'(A,1X,I4,A,I4,1X,I4)') 'bond', bonds_array(1,i), & + '; type = ', bonds_array(2,i), bonds_array(3,i) END DO END PROGRAM bonds @@ -1465,6 +1520,207 @@ Procedures Bound to the :f:type:`lammps` Derived Type -------- +.. f:subroutine:: gather(self, name, count, data) + + Gather the named per-atom, per-atom fix, per-atom compute, or fix + property/atom-based entities from all processes, in order by atom ID. + + .. versionadded:: TBD + + This subroutine gathers data from all processes and stores them in a + one-dimensional allocatable array. The array *data* will be + ordered by atom ID, which requires consecutive IDs (1 to *natoms*\ ). If you + need a similar array but for non-consecutive atom IDs, see + :cpp:func:`lammps_gather_concat`; for a similar array but for a subset of + atoms, see :cpp:func:`lammps_gather_subset`. + + The *data* array will be ordered in groups of *count* values, sorted by atom + ID (e.g., if *name* is *x*, then *data* is [x(1,1), x(2,1), x(3,1), x(1,2), + x(2,2), x(3,2), x(1,3), :math:`\dots`]); *data* must be ``ALLOCATABLE`` and + will be allocated to length (*count*\ :math:`{}\times{}`\ *natoms*), as + queried by :f:func:`get_natoms`. + + This function will return an error if fix or compute data are requested and + the fix or compute ID given does not have per-atom data. See the note about + re-interpreting the vector as a matrix at :f:subr:`gather_atoms`. + + This function is not compatible with ``-DLAMMPS_BIGBIG``. + + :p character(len=\*) name: desired quantity (e.g., "x" or "mask" for atom + properties, "f_id" for per-atom fix data, "c_id" for per-atom compute data, + "d_name" or "i_name" for fix property/atom vectors with *count* = 1, + "d2_name" or "i2_name" for fix property/atom vectors with + *count*\ :math:`{}> 1`) + :p integer(c_int) count: number of per-atom values (e.g., 1 for *type* or + *charge*, 3 for *x* or *f*); use *count* = 3 with *image* if you want the + image flags unpacked into (*x*,\ *y*,\ *z*) components. + :p real(c_double) data [dimension(:),allocatable]: array into which to store + the data. Array *must* have the ``ALLOCATABLE`` attribute and be of rank 1 + (i.e., ``DIMENSION(:)``). If this array is already allocated, it will be + reallocated to fit the length of the incoming data. + :to: :cpp:func:`lammps_gather` + +-------- + +.. f:subroutine:: gather_concat(self, name, count, data) + + Gather the named per-atom, per-atom fix, per-atom compute, or fix + property/atom-based entities from all processes, unordered. + + .. versionadded:: TBD + + This subroutine gathers data for all atoms and stores them in a + one-dimensional allocatable array. The data will be a + concatenation of chunks from each processor's owned atoms, in whatever order + the atoms are in on each processor. This process has no requirement that the + atom IDs be consecutive. If you need the ID of each atom, you can do another + call to either :f:subr:`gather_atoms_concat` or :f:subr:`gather_concat` with + *name* set to ``id``. If you have consecutive IDs and want the data to be in + order, use :f:subr:`gather`; for a similar array but for a subset of + atoms, use :f:subr:`gather_subset`. + + The *data* array will be in groups of *count* values, with *natoms* groups + total, but not in order by atom ID (e.g., if *name* is *x* and *count* is 3, + then *data* might be something like [x(1,11), x(2,11), x(3,11), x(1,3), + x(2,3), x(3,3), x(1,5), :math:`\dots`]); *data* must be ``ALLOCATABLE`` and + will be allocated to length (*count* :math:`\times` *natoms*), as queried by + :f:func:`get_natoms`. + + This function is not compatible with ``-DLAMMPS_BIGBIG``. + + :p character(len=\*) name: desired quantity (e.g., "x" or "mask" for atom + properties, "f_id" for per-atom fix data, "c_id" for per-atom compute data, + "d_name" or "i_name" for fix property/atom vectors with *count* = 1, + "d2_name" or "i2_name" for fix property/atom vectors with + *count*\ :math:`{}> 1`) + :p integer(c_int) count: number of per-atom values you expect per atom + (e.g., 1 for *type*, *mask*, or *charge*; 3 for *x*, *v*, or *f*). Use + *count* = 3 with *image* if you want a single image flag unpacked into + *x*/*y*/*z* components. + :p data: array into which to store + the data. Array *must* have the ``ALLOCATABLE`` attribute and be of rank 1 + (i.e., ``DIMENSION(:)``). If this array is already allocated, it will be + reallocated to fit the length of the incoming data. It should have type + ``INTEGER(c_int)`` if expecting integer data and ``REAL(c_double)`` if + expecting floating-point data. + :ptype data: polymorphic,dimension(:),allocatable + :to: :cpp:func:`lammps_gather_concat` + +-------- + +.. f:subroutine:: gather_subset(name, count, ids, data) + + Gather the named per-atom, per-atom fix, per-atom compute, or fix + property/atom-based entities from all processes for a subset of atoms. + + .. versionadded:: TBD + + This subroutine gathers data for the requested atom IDs and stores them in a + one-dimensional allocatable array. The data will be ordered by atom ID, but + there is no requirement that the IDs be consecutive. If you wish to return a + similar array for *all* the atoms, use :f:subr:`gather` or + :f:subr:`gather_concat`. + + The *data* array will be in groups of *count* values, sorted by atom ID in + the same order as the array *ids* (e.g., if *name* is *x*, *count* = 3, and + *ids* is [100, 57, 210], then *data* might look like [*x*\ (1,100), + *x*\ (2,100), *x*\ (3,100), *x*\ (1,57), *x*\ (2,57), *x*\ (3,57), + *x*\ (1,210), :math:`\dots`]); *ids* must be provided by the user, and + *data* must have the ``ALLOCATABLE`` attribute and be of rank 1 (i.e., + ``DIMENSION(:)``). If *data* is already allocated, it will be reallocated to + fit the length of the incoming data. + + This function is not compatible with ``-DLAMMPS_BIGBIG``. + + :p character(len=\*) name: quantity to be scattered + + :p integer(c_int) ids [dimension(:)]: atom IDs corresponding to the atoms + being scattered (e.g., "x" or "f" for atom properties, "f_id" for per-atom + fix data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix + property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix + property/atom vectors with *count*\ :math:`{} > 1`) + :p integer(c_int) count: number of per-atom values you expect per atom + (e.g., 1 for *type*, *mask*, or *charge*; 3 for *x*, *v*, or *f*). Use + *count* = 3 with *image* if you want a single image flag unpacked into + *x*/*y*/*z* components. + :p data: per-atom values packed into a one-dimensional array containing the + data to be scattered. This array must have the ``ALLOCATABLE`` attribute + and will be allocated either to the same length as *ids* + (for *mask*, *type*, etc.) or to three times its length (for *x*, *f*, + etc.); the array must be rank 1 and be of type ``INTEGER(c_int)`` (e.g., + for *mask* or *type*) or of type ``REAL(c_double)`` (e.g., for *charge*, + *x*, or *f*). + :ptype data: polymorphic,dimension(:),allocatable + :to: :cpp:func:`lammps_gather_subset` + +-------- + +.. f:subroutine:: scatter(name, data) + + This function calls :cpp:func:`lammps_scatter` to scatter the named + per-atom, per-atom fix, per-atom compute, or fix property/atom-based entity + in *data* to all processes. + + .. versionadded:: TBD + + This subroutine takes data stored in a one-dimensional array supplied by the + user and scatters them to all atoms on all processes. The data must be + ordered by atom ID, with the requirement that the IDs be consecutive. Use + :f:subr:`scatter_subset` to scatter data for some (or all) atoms, unordered. + + The *data* array needs to be ordered in groups of *count* values, sorted by + atom ID (e.g., if *name* is *x* and *count* = 3, then *data* = [*x*\ (1,1), + *x*\ (2,1), *x*\ (3,1), *x*\ (1,2), *x*\ (2,2), *x*\ (3,2), *x*\ (1,3), + :math:`\dots`]); *data* must be of length (*count* :math:`\times` *natoms*). + + This function is not compatible with ``-DLAMMPS_BIGBIG``. + + :p character(len=*) name: desired quantity (e.g., "x" or "f" for atom + properties, "f_id" for per-atom fix data, "c_id" for per-atom compute data, + "d_name" or "i_name" for fix property/atom vectors with *count* = 1, + "d2_name" or "i2_name" for fix property/atom vectors with *count*\ + :math:`{} > 1`) + :p data: per-atom values packed in a one-dimensional array; *data* should be + of type ``INTEGER(c_int)`` or ``REAL(c_double)``, depending on the type of + data being scattered, and be of rank 1 (i.e., ``DIMENSION(:)``). + :ptype data: polymorphic,dimension(:) + :to: :cpp:func:`lammps_scatter` + +-------- + +.. f:subroutine:: scatter_subset(name, ids, data) + + This function calls :cpp:func:`lammps_scatter_subset` to scatter the named + per-atom, per-atom fix, per-atom compute, or fix property/atom-based + entities in *data* from a subset of atoms to all processes. + + .. versionadded:: TBD + + This subroutine takes data stored in a one-dimensional array supplied by the + user and scatters them to a subset of atoms on all processes. The array + *data* contains data associated with atom IDs, but there is no requirement + that the IDs be consecutive, as they are provided in a separate array. + Use :f:subr:`scatter` to scatter data for all atoms, in order. + + The *data* array needs to be organized in groups of *count* values, with the + groups in the same order as the array *ids*. For example, if you want *data* + to be the array [x(1,1), x(2,1), x(3,1), x(1,100), x(2,100), x(3,100), + x(1,57), x(2,57), x(3,57)], then *count* = 3 and *ids* = [1, 100, 57]. + + This function is not compatible with ``-DLAMMPS_BIGBIG``. + + :p character(len=\*) name: desired quantity (e.g., "x" or "mask" for atom + properties, "f_id" for per-atom fix data, "c_id" for per-atom compute data, + "d_name" or "i_name" for fix property/atom vectors with *count* = 1, + "d2_name" or "i2_name" for fix property/atom vectors with + *count*\ :math:`{}> 1`) + :p integer(c_int) ids: list of atom IDs to scatter data for + :p polymorphic data [dimension(:)]: per-atom values packed in a + one-dimensional array of length *size(ids)* \* *count*. + :to: :cpp:func:`lammps_scatter_subset` + +-------- + .. f:subroutine:: create_atoms([id,] type, x, [v,] [image,] [bexpand]) This method calls :cpp:func:`lammps_create_atoms` to create additional atoms @@ -1475,24 +1731,25 @@ Procedures Bound to the :f:type:`lammps` Derived Type :p integer(c_int) type [dimension(N)]: vector of :math:`N` atom types (required/see note below) - :p real(c_double) x [dimension(3N)]: vector of :math:`3N` x/y/z positions + :p real(c_double) x [dimension(3N)]: vector of :math:`3N\ x/y/z` positions of the new atoms, arranged as :math:`[x_1,y_1,z_1,x_2,y_2,\dotsc]` (required/see note below) - :o integer(kind=\*) id [dimension(N)]: vector of :math:`N` atom IDs; if - absent, LAMMPS will generate them for you. \*The ``KIND`` parameter should - be ``c_int`` unless LAMMPS was compiled with ``-DLAMMPS_BIGBIG``, in which - case it should be ``c_int64_t``. - :o real(c_double) v [dimension(3N)]: vector of :math:`3N` x/y/z velocities - of the new atoms, arranged as :math:`[v_{1,x},v_{1,y},v_{1,z},v_{2,x}, - \dotsc]`; if absent, they will be set to zero - :o integer(kind=\*) image [dimension(N)]: vector of :math:`N` image flags; - if absent, they are set to zero. \*The ``KIND`` parameter should be + :o integer(kind=\*) id [dimension(N),optional]: vector of :math:`N` atom + IDs; if absent, LAMMPS will generate them for you. \*The ``KIND`` parameter + should be ``c_int`` unless LAMMPS was compiled with ``-DLAMMPS_BIGBIG``, in + which case it should be ``c_int64_t``. + :o real(c_double) v [dimension(3N),optional]: vector of :math:`3N` + *x*\ /*y*\ /*z* velocities of the new atoms, arranged as + :math:`[v_{1,x},v_{1,y},v_{1,z},v_{2,x}, \dotsc]`; if absent, they will be + set to zero + :o integer(kind=\*) image [dimension(N),optional]: vector of :math:`N` image + flags; if absent, they are set to zero. \*The ``KIND`` parameter should be ``c_int`` unless LAMMPS was compiled with ``-DLAMMPS_BIGBIG``, in which case it should be ``c_int64_t``. See note below. :o bexpand: if ``.TRUE.``, atoms outside of shrink-wrap boundaries will be created, not dropped, and the box dimensions will be extended. Default is ``.FALSE.`` - :otype bexpand: logical + :otype bexpand: logical,optional :to: :cpp:func:`lammps_create_atoms` .. note:: @@ -1541,11 +1798,11 @@ Procedures Bound to the :f:type:`lammps` Derived Type :p character(len=\*) style: String used to search for pair style instance. :o exact: Flag to control whether style should match exactly or only a regular expression/sub-string match is applied. Default: ``.TRUE.``. - :otype exact: logical - :o integer(c_int) nsub: Match *nsub*\ th hybrid sub-style instance of - the same style. Default: 0. - :o integer(c_int) reqid: Request ID to identify the neighbor list in - case there are multiple requests from the same pair style instance. + :otype exact: logical,optional + :o integer(c_int) nsub [optional]: Match *nsub*\ th hybrid sub-style + instance of the same style. Default: 0. + :o integer(c_int) reqid [optional]: Request ID to identify the neighbor list + in case there are multiple requests from the same pair style instance. Default: 0. :to: :cpp:func:`lammps_find_pair_neighlist` :r integer(c_int) index: Neighbor list index if found, otherwise @@ -1564,8 +1821,9 @@ Procedures Bound to the :f:type:`lammps` Derived Type fixes with multiple neighbor list requests. :p character(len=\*) id: Identifier of fix instance - :o integer(c_int) reqid: request ID to identify the neighbor list in cases - in which there are multiple requests from the same fix. Default: 0. + :o integer(c_int) reqid [optional]: request ID to identify the neighbor list + in cases in which there are multiple requests from the same fix. + Default: 0. :to: :cpp:func:`lammps_find_fix_neighlist` :r index: neighbor list index if found, otherwise :math:`-1` :rtype index: integer(c_int) @@ -1583,8 +1841,9 @@ Procedures Bound to the :f:type:`lammps` Derived Type in case a compute has multiple neighbor list requests. :p character(len=\*) id: Identifier of compute instance. - :o integer(c_int) reqid: request ID to identify the neighbor list in cases - in which there are multiple requests from the same compute. Default: 0. + :o integer(c_int) reqid [optional]: request ID to identify the neighbor list + in cases in which there are multiple requests from the same compute. + Default: 0. :to: :cpp:func:`lammps_find_compute_neighlist` :r index: neighbor list index if found, otherwise :math:`-1`. :rtype index: integer(c_int) @@ -1602,6 +1861,25 @@ Procedures Bound to the :f:type:`lammps` Derived Type :r inum: number of entries in neighbor list, or :math:`-1` if *idx* is not a valid index. :rtype inum: integer(c_int) + :to: :cpp:func:`lammps_neighlist_num_elements` + +-------- + +.. f:subroutine:: neighlist_element_neighbors(idx, element, iatom, neighbors) + + Return atom local index, number of neighbors, and array of neighbor local + atom indices of a neighbor list entry. + + .. versionadded:: 3Nov2022 + + :p integer(c_int) idx: index of this neighbor list in the list of all + neighbor lists + :p integer(c_int) element: index of this neighbor list entry + :p integer(c_int) iatom: local atom index (i.e., in the range + [1,nlocal+nghost]; -1 if invalid or element value + :p integer(c_int) neighbors [dimension(:),pointer]: pointer to an array of + neighboring atom local indices + :to: :cpp:func:`lammps_neighlist_element_neighbors` -------- @@ -1806,9 +2084,10 @@ Procedures Bound to the :f:type:`lammps` Derived Type .. versionadded:: 3Nov2022 - This function is analogous to the :py:func`installed_packages` function in - the Python API. The optional argument *length* sets the length of each - string in the vector *package* (default: 31). + This function is analogous to the :py:meth:`installed_packages + ` function in the Python API. + The optional argument *length* sets the length of each string in the vector + *package* (default: 31). :p character(len=:) package [dimension(:),allocatable]: list of packages; *must* have the ``ALLOCATABLE`` attribute and be of rank 1 @@ -2044,11 +2323,12 @@ Procedures Bound to the :f:type:`lammps` Derived Type necessary to combine the values of three integers representing the image flags in the :math:`x`-, :math:`y`-, and :math:`z`-directions. Unless LAMMPS is compiled with ``-DLAMMPS_BIGBIG``, those integers are limited to 10-bit - signed integers :math:`[-512,512]`. If ``-DLAMMPS_BIGBIG`` was used when + signed integers :math:`[-512,512)`. If ``-DLAMMPS_BIGBIG`` was used when compiling, then the return value is of kind ``c_int64_t`` instead of kind ``c_int``, and the valid range for the individual image flags becomes - :math:`[-1048576,1048575]` (i.e., the range of a 21-bit signed integer). - There is no check on whether the arguments conform to these requirements. + :math:`[-1048576,1048575)` (i.e., the range of a 21-bit signed integer). + There is no check on whether the arguments conform to these requirements; + values out of range will simply be wrapped back into the interval. :p integer(c_int) ix: image flag in :math:`x`-direction :p integer(c_int) iy: image flag in :math:`y`-direction @@ -2096,6 +2376,355 @@ Procedures Bound to the :f:type:`lammps` Derived Type -------- +.. f:subroutine:: set_fix_external_callback(id, callback, caller) + + Set the callback function for a :doc:`fix external ` instance + with the given ID. + + .. versionadded:: TBD + + Fix :doc:`external ` allows programs that are running LAMMPS + through its library interface to modify certain LAMMPS properties on + specific time steps, similar to the way other fixes do. + + This subroutine sets the callback function for use with the "pf/callback" + mode. The function should have Fortran language bindings with the following + interface, which depends on how LAMMPS was compiled: + + .. code-block:: Fortran + + ABSTRACT INTERFACE + SUBROUTINE external_callback(caller, timestep, ids, x, fexternal) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double, c_int64_t + CLASS(*), INTENT(INOUT) :: caller + INTEGER(c_bigint), INTENT(IN) :: timestep + INTEGER(c_tagint), DIMENSION(:), INTENT(IN) :: ids + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: fexternal + END SUBROUTINE external_callback + END INTERFACE + + where ``c_bigint`` is ``c_int`` if ``-DLAMMPS_SMALLSMALL`` was used and + ``c_int64_t`` otherwise; and ``c_tagint`` is ``c_int64_t`` if + ``-DLAMMPS_BIGBIG`` was used and ``c_int`` otherwise. + + The argument *caller* to :f:subr:`set_fix_external_callback` is unlimited + polymorphic (i.e., it can be any Fortran object you want to pass to the + calling function) and will be available as the first argument to the + callback function. It can be your LAMMPS instance, which you might need if + the callback function needs access to the library interface. + The argument must be a scalar; to pass non-scalar data, wrap those data in + a derived type and pass an instance of the derived type to *caller*. + + The array *ids* is an array of length *nlocal* (as accessed from the + :cpp:class:`Atom` class or through :f:func:`extract_global`). The arrays + *x* and *fexternal* are :math:`3 \times {}`\ *nlocal* arrays; these are + transposed from what they would look like in C (see note about array index + order at :f:func:`extract_atom`). + + The callback mechanism is one of two ways that forces can be applied to a + simulation with the help of :doc:`fix external `. The + alternative is *array* mode, where one calls + :f:func:`fix_external_get_force`. + + Please see the documentation for :doc:`fix external ` for + more information about how to use the fix and couple it with external + programs. + + :p character(len=*) id: ID of :doc:`fix external ` instance + :p callback: subroutine :doc:`fix external ` should call + :ptype callback: external + :p class(*) caller [optional]: object you wish to pass to the callback + procedure (must be a scalar; see note) + :to: :cpp:func:`lammps_set_fix_external_callback` + + .. note:: + + The interface for your callback function must match types precisely + with the abstract interface block given above. **The compiler probably + will not be able to check this for you.** In particular, the first + argument ("caller") must be of type ``CLASS(*)`` or you will probably + get a segmentation fault or at least a misinterpretation of whatever is + in memory there. You can resolve the object using the ``SELECT TYPE`` + construct. An example callback function (assuming LAMMPS was compiled + with ``-DLAMMPS_SMALLBIG``) that applies something akin to Hooke's Law + (with each atom having a different *k* value) is shown below. + + .. code-block:: Fortran + + MODULE stuff + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double, c_int64_t + USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY : error_unit + IMPLICIT NONE + + TYPE shield + REAL(c_double), DIMENSION(:,:), ALLOCATABLE :: k + ! assume k gets allocated to dimension(3,nlocal) at some point + ! and assigned values + END TYPE shield + + SUBROUTINE my_callback(caller, timestep, ids, x, fexternal) + CLASS(*), INTENT(INOUT) :: caller + INTEGER(c_int), INTENT(IN) :: timestep + INTEGER(c_int64_t), INTENT(IN) :: ids + REAL(c_double), INTENT(IN) :: x(:,:) + REAL(c_double), INTENT(OUT) :: fexternal(:,:) + + SELECT TYPE (caller) + TYPE IS (shield) + fexternal = - caller%k * x + CLASS DEFAULT + WRITE(error_unit,*) 'UH OH...' + END SELECT + END SUBROUTINE my_callback + END MODULE stuff + + ! then, when assigning the callback function, do this: + PROGRAM example + USE LIBLAMMPS + USE stuff + TYPE(lammps) :: lmp + TYPE(shield) :: my_shield + lmp = lammps() + CALL lmp%command('fix ext all external pf/callback 1 1') + CALL lmp%set_fix_external_callback('ext', my_callback, my_shield) + END PROGRAM example + +-------- + +.. f:function:: fix_external_get_force(id) + + Get pointer to the force array storage in a fix external instance with the + given ID. + + .. versionadded:: TBD + + Fix :doc:`external ` allows programs that are running LAMMPS + through its library interfaces to add or modify certain LAMMPS properties on + specific time steps, similar to the way other fixes do. + + This function provides access to the per-atom force storage in a fix + external instance with the given fix-ID to be added to the individual + atoms when using the "pf/array" mode. The *fexternal* array can be + accessed like other "native" per-atom arrays accessible via the + :f:func:`extract_atom` function. Please note that the array + stores the forces for *local* atoms for each MPI rank, in the order + determined by the neighbor list build. Because the underlying data + structures can change as well as the order of atom as they migrate between + MPI processes because of the domain decomposition parallelization, this + function should be always called immediately before the forces are going to + be set to get an up-to-date pointer. You can use, for example, + :f:func:`extract_setting` to obtain the number of local atoms `nlocal` and + then assume the dimensions of the returned force array as + ``REAL(c_double) :: force(3,nlocal)``. + + This function is an alternative to the callback mechanism in fix external + set up by :f:subr:`set_fix_external_callback`. The main difference is that + this mechanism can be used when forces are to be pre-computed and the + control alternates between LAMMPS and the external driver, while the + callback mechanism can call an external subroutine to compute the force when + the fix is triggered and needs them. + + Please see the documentation for :doc:`fix external ` for + more information about how to use the fix and how to couple it with an + external program. + + :p character(len=*) id: ID of :doc:`fix external ` instance + :to: :cpp:func:`lammps_fix_external_get_force` + :r fexternal: pointer to the per-atom force array allocated by the fix + :rtype fexternal: real(c_double),dimension(3,nlocal) + +-------- + +.. f:subroutine:: fix_external_set_energy_global(id, eng) + + Set the global energy contribution for a :doc:`fix external ` + instance with the given ID. + + .. versionadded:: TBD + + This is a companion function to :f:func:`set_fix_external_callback` + and :f:func:`fix_external_get_force` that also sets the contribution to the + global energy from the external program. The value of the *eng* argument + will be stored in the fix and applied on the current and all following + time steps until changed by another call to this function. The energy is in + energy units as determined by the current :doc:`units ` settings and + is the **total** energy of the contribution. Thus, when running in + parallel, all MPI processes have to call this function with the **same** + value, and this will be returned as a scalar property of the fix external + instance when accessed in LAMMPS input commands or from variables. + + Please see the documentation for :doc:`fix external ` for more + information about how to use the fix and how to couple it with an external + program. + + :p character(len=*) id: fix ID of fix external instance + :p real(c_double) eng: total energy to be added to the global energy + :to: :cpp:func:`lammps_fix_external_set_energy_global` + +-------- + +.. f:subroutine:: fix_external_set_virial_global(id, virial) + + Set the global virial contribution for a fix external instance with the + given ID. + + .. versionadded:: TBD + + This is a companion function to :f:subr:`set_fix_external_callback` + and :f:func:`fix_external_get_force` to set the contribution to the global + virial from an external program. + + The six values of the *virial* array will be stored in the fix and applied + on the current and all following time steps until changed by another call + to this function. The components of the virial need to be stored in the + following order: *xx*, *yy*, *zz*, *xy*, *xz*, *yz*. In LAMMPS, the virial + is stored internally as `stress*volume` in units of `pressure*volume` as + determined by the current :doc:`units ` settings and is the + **total** contribution. Thus, when running in parallel, all MPI processes + have to call this function with the **same** value, and this will then be + added by fix external. + + Please see the documentation for :doc:`fix external ` for + more information about how to use the fix and how to couple it with an + external code. + + :p character(len=*) id: fix ID of fix external instance + :p real(c_double) virial [dimension(6)]: the six global stress tensor + components to be added to the global virial + :to: :cpp:func:`lammps_fix_external_set_virial_global` + +-------- + +.. f:subroutine:: fix_external_set_energy_peratom(id, eng) + + Set the per-atom energy contribution for a fix external instance with the + given ID. + + .. versionadded:: TBD + + This is a companion function to :f:subr:`set_fix_external_callback` to set + the per-atom energy contribution due to the fix from the external program as + part of the callback function. For this to work, the LAMMPS object must be + passed as part of the *caller* argument when registering the callback + function, or the callback function must otherwise have access to the + LAMMPS object, such as through a module-based pointer. + + .. note:: + + This function is fully independent from + :f:subr:`fix_external_set_energy_global` and will **NOT** add any + contributions to the global energy tally and will **NOT** check whether + the sum of the contributions added here are consistent with the global + added energy. + + Please see the documentation for :doc:`fix external ` for + more information about how to use the fix and how to couple it with an + external code. + + :p character(len=*) id: fix ID of the fix external instance + :p real(c_double) eng [dimension(:)]: array of length *nlocal* containing + the energy to add to the per-atom energy + :to: :cpp:func:`lammps_fix_external_set_energy_peratom` + +-------- + +.. f:subroutine:: set_fix_external_set_virial_peratom(id, virial) + + This is a companion function to :f:subr:`set_fix_external_callback` to set + the per-atom virial contribution due to the fix from the external program as + part of the callback function. For this to work, the LAMMPS object must be + passed as the *caller* argument when registering the callback function. + + .. versionadded:: TBD + + .. note:: + + This function is fully independent from + :f:subr:`fix_external_set_virial_global` and will **NOT** add any + contributions to the global virial tally and **NOT** check whether the + sum of the contributions added here are consistent with the global added + virial. + + The order and units of the per-atom stress tensor elements are the same + as for the global virial. The type and dimensions of the per-atom virial + array must be ``REAL(c_double), DIMENSION(6,nlocal)``. + + Please see the documentation for :doc:`fix external ` for + more information about how to use the fix and how to couple it with an + external program. + + :p character(len=*) id: fix ID of fix external instance + :p real(c_double) virial [dimension(:,:)]: an array of :math:`6 \times{}`\ + *nlocal* components to be added to the per-atom virial + :to: :cpp:func:`lammps_set_virial_peratom` + +-------- + +.. f:subroutine:: fix_external_set_vector_length(id, length) + + Set the vector length for a global vector stored with fix external for + analysis. + + .. versionadded:: TBD + + This is a companion function to :f:subr:`set_fix_external_callback` and + :f:func:`fix_external_get_force` to set the length of a global vector of + properties that will be stored with the fix via + :f:subr:`fix_external_set_vector`. + + This function needs to be called **before** a call to + :f:subr:`fix_external_set_vector` and **before** a run or minimize command. + When running in parallel, it must be called from **all** MPI + processes with the same length argument. + + Please see the documentation for :doc:`fix external ` for + more information about how to use the fix and how to couple it with an + external program. + + :p character(len=*) id: fix ID of fix external instance + :p integer(c_int) length: length of the global vector to be stored with the + fix + :to: :cpp:func:`lammps_fix_external_set_vector_length` + +-------- + +.. f:subroutine:: fix_external_set_vector(id, idx, val) + + Store a global vector value for a fix external instance with the given ID. + + .. versionadded:: TBD + + This is a companion function to :f:subr:`set_fix_external_callback` and + :f:func:`fix_external_get_force` to set the values of a global vector of + properties that will be stored with the fix and can be accessed from + within LAMMPS input commands (e.g., fix ave/time or variables) when used + in a vector context. + + This function needs to be called **after** a call to + :f:subr:`fix_external_set_vector_length` and **before** a run or minimize + command. When running in parallel, it must be called from **all** MPI + processes with the **same** *idx* and *val* parameters. The variable + *val* is assumed to be extensive. + + .. note:: + + The index in the *idx* parameter is 1-based (i.e., the first element + is set with *idx*\ :math:`{} = 1`, and the last element of the vector + with *idx*\ :math:`{} = N`, where :math:`N` is the value of the *length* + parameter of the call to :f:subr:`fix_external_set_vector_length`). + + Please see the documentation for :doc:`fix external ` for + more information about how to use the fix and how to couple it with an + external code. + + :p character(len=*) id: ID of fix external instance + :p integer(c_int) idx: 1-based index in global vector + :p integer(c_int) val: value to be stored in global vector at index *idx* + :to: :cpp:func:`lammps_fix_external_set_vector` + +-------- + .. f:subroutine:: flush_buffers() This function calls :cpp:func:`lammps_flush_buffers`, which flushes buffered @@ -2130,7 +2759,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type .. versionadded:: 3Nov2022 This function can be used from signal handlers or multi-threaded - applications to cleanly terminate an ongoing run. + applications to terminate an ongoing run cleanly. :to: :cpp:func:`lammps_force_timeout` @@ -2169,12 +2798,11 @@ Procedures Bound to the :f:type:`lammps` Derived Type in the event of an error inside of LAMMPS that resulted in a :ref:`C++ exception `. A suitable buffer for a string has to be provided. If the internally-stored error message is longer than the - string and the string does not have ``ALLOCATABLE`` length, it will be - truncated accordingly. The optional argument *status* indicates the - kind of error: a "1" indicates an error that occurred on all MPI ranks and - is often recoverable, while a "2" indicates an abort that would happen only - in a single MPI rank and thus may not be recoverable, as other MPI ranks may - be waiting on the failing MPI rank(s) to send messages. + string, it will be truncated accordingly. The optional argument *status* + indicates the kind of error: a "1" indicates an error that occurred on all + MPI ranks and is often recoverable, while a "2" indicates an abort that + would happen only in a single MPI rank and thus may not be recoverable, as + other MPI ranks may be waiting on the failing MPI rank(s) to send messages. .. note:: diff --git a/doc/src/Library_utility.rst b/doc/src/Library_utility.rst index 9d16daae0b..e555b79c0b 100644 --- a/doc/src/Library_utility.rst +++ b/doc/src/Library_utility.rst @@ -7,6 +7,7 @@ functions. They do not directly call the LAMMPS library. - :cpp:func:`lammps_encode_image_flags` - :cpp:func:`lammps_decode_image_flags` - :cpp:func:`lammps_set_fix_external_callback` +- :cpp:func:`lammps_fix_external_get_force` - :cpp:func:`lammps_fix_external_set_energy_global` - :cpp:func:`lammps_fix_external_set_energy_peratom` - :cpp:func:`lammps_fix_external_set_virial_global` @@ -44,6 +45,11 @@ where such memory buffers were allocated that require the use of ----------------------- +.. doxygenfunction:: lammps_fix_external_get_force + :project: progguide + +----------------------- + .. doxygenfunction:: lammps_fix_external_set_energy_global :project: progguide diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index 3ab7a26d25..e25b5ec206 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -35,7 +35,7 @@ MODULE LIBLAMMPS USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_ptr, c_null_ptr, C_ASSOCIATED, & C_LOC, c_int, c_int64_t, c_char, c_null_char, c_double, c_size_t, & - C_F_POINTER + C_F_POINTER, c_funptr, C_FUNLOC IMPLICIT NONE PRIVATE @@ -55,7 +55,7 @@ MODULE LIBLAMMPS LAMMPS_DOUBLE_2D = 3, & ! two-dimensional 64-bit double array LAMMPS_INT64 = 4, & ! 64-bit integer (or array) LAMMPS_INT64_2D = 5, & ! two-dimensional 64-bit integer array - LAMMPS_STRING = 6, & ! C-String + LAMMPS_STRING = 6, & ! string LMP_STYLE_GLOBAL = 0, & ! request global compute/fix/etc. data LMP_STYLE_ATOM = 1, & ! request per-atom compute/fix/etc. data LMP_STYLE_LOCAL = 2, & ! request local compute/fix/etc. data @@ -64,7 +64,7 @@ MODULE LIBLAMMPS LMP_TYPE_ARRAY = 2, & ! request array LMP_SIZE_VECTOR = 3, & ! request size of vector LMP_SIZE_ROWS = 4, & ! request rows (actually columns) - LMP_SIZE_COLS = 5, & ! request colums (actually rows) + LMP_SIZE_COLS = 5, & ! request columns (actually rows) LMP_ERROR_WARNING = 0, & ! call Error::warning() LMP_ERROR_ONE = 1, & ! call Error::one() (from this MPI rank) LMP_ERROR_ALL = 2, & ! call Error::all() (from all MPI ranks) @@ -75,6 +75,9 @@ MODULE LIBLAMMPS LMP_VAR_VECTOR = 2, & ! vector variables LMP_VAR_STRING = 3 ! string variables (everything else) + ! Constants we set once (in the constructor) and never need to check again + INTEGER(c_int), SAVE :: SIZE_TAGINT, SIZE_BIGINT, SIZE_IMAGEINT + ! "Constants" to use with extract_compute and friends TYPE lammps_style INTEGER(c_int) :: global, atom, local @@ -132,7 +135,24 @@ MODULE LIBLAMMPS PROCEDURE, PRIVATE :: lmp_gather_bonds_big GENERIC :: gather_bonds => lmp_gather_bonds_small, & lmp_gather_bonds_big -! + PROCEDURE, PRIVATE :: lmp_gather_int + PROCEDURE, PRIVATE :: lmp_gather_double + GENERIC :: gather => lmp_gather_int, lmp_gather_double + PROCEDURE, PRIVATE :: lmp_gather_concat_int + PROCEDURE, PRIVATE :: lmp_gather_concat_double + GENERIC :: gather_concat => lmp_gather_concat_int, & + lmp_gather_concat_double + PROCEDURE, PRIVATE :: lmp_gather_subset_int + PROCEDURE, PRIVATE :: lmp_gather_subset_double + GENERIC :: gather_subset => lmp_gather_subset_int, & + lmp_gather_subset_double + PROCEDURE, PRIVATE :: lmp_scatter_int + PROCEDURE, PRIVATE :: lmp_scatter_double + GENERIC :: scatter => lmp_scatter_int, lmp_scatter_double + PROCEDURE, PRIVATE :: lmp_scatter_subset_int + PROCEDURE, PRIVATE :: lmp_scatter_subset_double + GENERIC :: scatter_subset => lmp_scatter_subset_int, & + lmp_scatter_subset_double PROCEDURE, PRIVATE :: lmp_create_atoms_int PROCEDURE, PRIVATE :: lmp_create_atoms_bigbig GENERIC :: create_atoms => lmp_create_atoms_int, & @@ -171,7 +191,19 @@ MODULE LIBLAMMPS PROCEDURE, PRIVATE :: lmp_decode_image_flags_bigbig GENERIC :: decode_image_flags => lmp_decode_image_flags, & lmp_decode_image_flags_bigbig -! + PROCEDURE :: set_fix_external_callback => lmp_set_fix_external_callback + PROCEDURE :: fix_external_get_force => lmp_fix_external_get_force + PROCEDURE :: fix_external_set_energy_global & + => lmp_fix_external_set_energy_global + PROCEDURE :: fix_external_set_virial_global & + => lmp_fix_external_set_virial_global + PROCEDURE :: fix_external_set_energy_peratom & + => lmp_fix_external_set_energy_peratom + PROCEDURE :: fix_external_set_virial_peratom & + => lmp_fix_external_set_virial_peratom + PROCEDURE :: fix_external_set_vector_length & + => lmp_fix_external_set_vector_length + PROCEDURE :: fix_external_set_vector => lmp_fix_external_set_vector PROCEDURE :: flush_buffers => lmp_flush_buffers PROCEDURE :: is_running => lmp_is_running PROCEDURE :: force_timeout => lmp_force_timeout @@ -261,6 +293,50 @@ MODULE LIBLAMMPS assign_int64_to_lammps_image_data END INTERFACE + ! Interface templates for fix external callbacks + ABSTRACT INTERFACE + SUBROUTINE external_callback_smallsmall(caller, timestep, ids, x, fexternal) + IMPORT :: c_int, c_double + CLASS(*), INTENT(INOUT) :: caller + INTEGER(c_int), INTENT(IN) :: timestep + INTEGER(c_int), DIMENSION(:), INTENT(IN) :: ids + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: fexternal + END SUBROUTINE external_callback_smallsmall + SUBROUTINE external_callback_smallbig(caller, timestep, ids, x, fexternal) + IMPORT :: c_int, c_double, c_int64_t + CLASS(*), INTENT(INOUT) :: caller + INTEGER(c_int64_t), INTENT(IN) :: timestep + INTEGER(c_int), DIMENSION(:), INTENT(IN) :: ids + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: fexternal + END SUBROUTINE external_callback_smallbig + SUBROUTINE external_callback_bigbig(caller, timestep, ids, x, fexternal) + IMPORT :: c_double, c_int64_t + CLASS(*), INTENT(INOUT) :: caller + INTEGER(c_int64_t), INTENT(IN) :: timestep + INTEGER(c_int64_t), DIMENSION(:), INTENT(IN) :: ids + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: fexternal + END SUBROUTINE external_callback_bigbig + END INTERFACE + + ! Derived type for fix external callback data + TYPE fix_external_data + CHARACTER(LEN=:), ALLOCATABLE :: id + PROCEDURE(external_callback_smallsmall), NOPASS, POINTER :: & + callback_smallsmall => NULL() + PROCEDURE(external_callback_smallbig), NOPASS, POINTER :: & + callback_smallbig => NULL() + PROCEDURE(external_callback_bigbig), NOPASS, POINTER :: & + callback_bigbig => NULL() + CLASS(*), POINTER :: caller => NULL() + CLASS(lammps), POINTER :: lammps_instance => NULL() + END TYPE fix_external_data + + ! Array used to store Fortran-facing callback functions for fix external + TYPE(fix_external_data), DIMENSION(:), ALLOCATABLE, TARGET, SAVE :: ext_data + ! interface definitions for calling functions in library.cpp INTERFACE FUNCTION lammps_open(argc, argv, comm) BIND(C,name='lammps_open_fortran') @@ -493,13 +569,42 @@ MODULE LIBLAMMPS TYPE(c_ptr), VALUE :: handle, data END SUBROUTINE lammps_gather_bonds - !SUBROUTINE lammps_gather + SUBROUTINE lammps_gather(handle, name, type, count, data) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, name, data + INTEGER(c_int), VALUE :: type, count + END SUBROUTINE lammps_gather - !SUBROUTINE lammps_gather_concat + SUBROUTINE lammps_gather_concat(handle, name, type, count, data) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, name, data + INTEGER(c_int), VALUE :: type, count + END SUBROUTINE lammps_gather_concat - !SUBROUTINE lammps_gather_subset + SUBROUTINE lammps_gather_subset(handle, name, type, count, ndata, ids, & + data) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, name, ids, data + INTEGER(c_int), VALUE :: type, count, ndata + END SUBROUTINE lammps_gather_subset - !SUBROUTINE lammps_scatter_subset + SUBROUTINE lammps_scatter(handle, name, type, count, data) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, name, data + INTEGER(c_int), VALUE :: type, count + END SUBROUTINE lammps_scatter + + SUBROUTINE lammps_scatter_subset(handle, name, type, count, ndata, ids, & + data) BIND(C) + IMPORT :: c_ptr, c_int + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, name, ids, data + INTEGER(c_int), VALUE :: count, ndata, type + END SUBROUTINE lammps_scatter_subset FUNCTION lammps_create_atoms(handle, n, id, type, x, v, image, bexpand) & BIND(C) @@ -705,15 +810,55 @@ MODULE LIBLAMMPS ! It is re-written in Fortran below. It was easier to do the same for ! lammps_decode_image_flags's equivalent. - !SUBROUTINE lammps_set_fix_external_callback ! may have trouble.... - !FUNCTION lammps_fix_external_get_force() ! returns real(c_double)(:) + SUBROUTINE lammps_set_fix_external_callback(handle, id, funcptr, ptr) & + BIND(C) + IMPORT :: c_ptr, c_funptr + TYPE(c_ptr), VALUE :: handle, id, ptr + TYPE(c_funptr), VALUE :: funcptr + END SUBROUTINE lammps_set_fix_external_callback - !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 + FUNCTION lammps_fix_external_get_force(handle, id) BIND(C) + IMPORT :: c_ptr + TYPE(c_ptr), VALUE :: handle, id + TYPE(c_ptr) :: lammps_fix_external_get_force + END FUNCTION lammps_fix_external_get_force + + SUBROUTINE lammps_fix_external_set_energy_global(handle, id, eng) BIND(C) + IMPORT :: c_ptr, c_double + TYPE(c_ptr), VALUE :: handle, id + REAL(c_double), VALUE :: eng + END SUBROUTINE lammps_fix_external_set_energy_global + + SUBROUTINE lammps_fix_external_set_virial_global(handle, id, virial) & + BIND(C) + IMPORT :: c_ptr + TYPE(c_ptr), VALUE :: handle, id, virial + END SUBROUTINE lammps_fix_external_set_virial_global + + SUBROUTINE lammps_fix_external_set_energy_peratom(handle, id, eng) BIND(C) + IMPORT :: c_ptr + TYPE(c_ptr), VALUE :: handle, id, eng + END SUBROUTINE lammps_fix_external_set_energy_peratom + + SUBROUTINE lammps_fix_external_set_virial_peratom(handle, id, virial) & + BIND(C) + IMPORT :: c_ptr + TYPE(c_ptr), VALUE :: handle, id, virial + END SUBROUTINE lammps_fix_external_set_virial_peratom + + SUBROUTINE lammps_fix_external_set_vector_length(handle, id, length) & + BIND(C) + IMPORT :: c_ptr, c_int + TYPE(c_ptr), VALUE :: handle, id + INTEGER(c_int), VALUE :: length + END SUBROUTINE lammps_fix_external_set_vector_length + + SUBROUTINE lammps_fix_external_set_vector(handle, id, idx, val) BIND(C) + IMPORT :: c_ptr, c_int, c_double + TYPE(c_ptr), VALUE :: handle, id + INTEGER(c_int), VALUE :: idx + REAL(c_double), VALUE :: val + END SUBROUTINE lammps_fix_external_set_vector SUBROUTINE lammps_flush_buffers(handle) BIND(C) IMPORT :: c_ptr @@ -784,28 +929,28 @@ CONTAINS INTEGER(c_int) :: i, c_comm, argc IF (PRESENT(args)) THEN - ! convert fortran argument list to c style - argc = SIZE(args) - ALLOCATE(argv(argc)) - DO i=1, argc - argv(i) = f2c_string(args(i)) - END DO + ! convert fortran argument list to c style + argc = SIZE(args) + ALLOCATE(argv(argc)) + DO i=1, argc + argv(i) = f2c_string(args(i)) + END DO ELSE - argc = 1 - ALLOCATE(argv(1)) - argv(1) = f2c_string("liblammps") + argc = 1 + ALLOCATE(argv(1)) + argv(1) = f2c_string("liblammps") ENDIF IF (PRESENT(comm)) THEN - c_comm = comm - lmp_open%handle = lammps_open(argc, argv, c_comm) + c_comm = comm + lmp_open%handle = lammps_open(argc, argv, c_comm) ELSE - lmp_open%handle = lammps_open_no_mpi(argc, argv, c_null_ptr) + lmp_open%handle = lammps_open_no_mpi(argc, argv, c_null_ptr) END IF ! Clean up allocated memory DO i=1, argc - CALL lammps_free(argv(i)) + CALL lammps_free(argv(i)) END DO DEALLOCATE(argv) @@ -816,6 +961,11 @@ CONTAINS lmp_open%type%scalar = LMP_TYPE_SCALAR lmp_open%type%vector = LMP_TYPE_VECTOR lmp_open%type%array = LMP_TYPE_ARRAY + + ! Assign constants for bigint and tagint for use elsewhere + SIZE_TAGINT = lmp_extract_setting(lmp_open, 'tagint') + SIZE_BIGINT = lmp_extract_setting(lmp_open, 'bigint') + SIZE_IMAGEINT = lmp_extract_setting(lmp_open, 'imageint') END FUNCTION lmp_open ! Combined Fortran wrapper around lammps_close() and lammps_mpi_finalize() @@ -826,10 +976,10 @@ CONTAINS CALL lammps_close(self%handle) IF (PRESENT(finalize)) THEN - IF (finalize) THEN - CALL lammps_kokkos_finalize() - CALL lammps_mpi_finalize() - END IF + IF (finalize) THEN + CALL lammps_kokkos_finalize() + CALL lammps_mpi_finalize() + END IF END IF END SUBROUTINE lmp_close @@ -997,7 +1147,7 @@ CONTAINS length = 3 CASE DEFAULT length = 1 - ! string cases doesn't use "length" + ! string cases do not use "length" END SELECT Cname = f2c_string(name) @@ -1518,7 +1668,7 @@ CONTAINS IF (count /= 1 .AND. count /= 3) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & 'gather_atoms_subset requires "count" to be 1 or 3 & - &[Fortran/gather_atoms]') + &[Fortran/gather_atoms_subset]') END IF ndata = SIZE(ids, KIND=c_int) @@ -1548,7 +1698,7 @@ CONTAINS IF (count /= 1 .AND. count /= 3) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & 'gather_atoms_subset requires "count" to be 1 or 3 & - &[Fortran/gather_atoms]') + &[Fortran/gather_atoms_subset]') END IF ndata = SIZE(ids, KIND=c_int) @@ -1642,7 +1792,8 @@ CONTAINS Ccount = SIZE(data, KIND=c_int) / Cndata IF (Ccount /= 1 .AND. Ccount /= 3) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & - 'scatter_atoms_subset requires either 1 or 3 data per atom') + 'scatter_atoms_subset requires either 1 or 3 data per atom & + &[Fortran/scatter_atoms_subset]') END IF Cname = f2c_string(name) @@ -1667,7 +1818,8 @@ CONTAINS Ccount = SIZE(data, KIND=c_int) / Cndata IF (Ccount /= 1 .AND. Ccount /= 3) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & - 'scatter_atoms_subset requires either 1 or 3 data per atom') + 'scatter_atoms_subset requires either 1 or 3 data per atom & + &[Fortran/scatter_atoms_subset]') END IF Cname = f2c_string(name) @@ -1682,19 +1834,16 @@ CONTAINS SUBROUTINE lmp_gather_bonds_small(self, data) CLASS(lammps), INTENT(IN) :: self INTEGER(c_int), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data - INTEGER(c_int) :: size_tagint, size_bigint INTEGER(c_int), POINTER :: nbonds_small INTEGER(c_int64_t), POINTER :: nbonds_big TYPE(c_ptr) :: Cdata - size_tagint = lmp_extract_setting(self, 'tagint') - IF (size_tagint /= 4_c_int) THEN + IF (SIZE_TAGINT /= 4_c_int) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & - 'Incompatible integer kind in gather_bonds [Fortran API]') + 'Incompatible integer kind in gather_bonds [Fortran/gather_bonds]') END IF IF (ALLOCATED(data)) DEALLOCATE(data) - size_bigint = lmp_extract_setting(self, 'bigint') - IF (size_bigint == 4_c_int) THEN + IF (SIZE_BIGINT == 4_c_int) THEN nbonds_small = lmp_extract_global(self, 'nbonds') ALLOCATE(data(3*nbonds_small)) ELSE @@ -1709,14 +1858,12 @@ CONTAINS SUBROUTINE lmp_gather_bonds_big(self, data) CLASS(lammps), INTENT(IN) :: self INTEGER(c_int64_t), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data - INTEGER(c_int) :: size_tagint INTEGER(c_int64_t), POINTER :: nbonds TYPE(c_ptr) :: Cdata - size_tagint = lmp_extract_setting(self, 'tagint') - IF (size_tagint /= 8_c_int) THEN + IF (SIZE_TAGINT /= 8_c_int) THEN CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & - 'Incompatible integer kind in gather_bonds [Fortran API]') + 'Incompatible integer kind in gather_bonds [Fortran/gather_bonds]') END IF nbonds = lmp_extract_global(self, 'nbonds') IF (ALLOCATED(data)) DEALLOCATE(data) @@ -1725,6 +1872,314 @@ CONTAINS CALL lammps_gather_bonds(self%handle, Cdata) END SUBROUTINE lmp_gather_bonds_big + ! equivalent function to lammps_gather (for int data) + SUBROUTINE lmp_gather_int(self, name, count, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), INTENT(IN) :: count + INTEGER(c_int), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data + TYPE(c_ptr) :: Cdata, Cname + INTEGER(c_int) :: natoms + INTEGER(c_int), PARAMETER :: Ctype = 0_c_int + REAL(c_double) :: dnatoms + CHARACTER(LEN=100) :: error_msg + + IF (count /= 1 .AND. count /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'gather requires "count" to be 1 or 3 [Fortran/gather]') + END IF + + dnatoms = lmp_get_natoms(self) + IF (dnatoms > HUGE(1_c_int)) THEN + WRITE(error_msg,'(A,1X,I0,1X,A)') & + 'Cannot use library function gather with more than', & + HUGE(0_c_int), 'atoms [Fortran/gather]' + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, error_msg) + END IF + natoms = NINT(dnatoms, c_int) + + Cname = f2c_string(name) + IF (ALLOCATED(data)) DEALLOCATE(data) + ALLOCATE(data(natoms*count)) + Cdata = C_LOC(data(1)) + CALL lammps_gather(self%handle, Cname, Ctype, count, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_gather_int + + ! equivalent function to lammps_gather_atoms (for doubles) + SUBROUTINE lmp_gather_double(self, name, count, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), INTENT(IN) :: count + REAL(c_double), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data + TYPE(c_ptr) :: Cdata, Cname + INTEGER(c_int) :: natoms + INTEGER(c_int), PARAMETER :: Ctype = 1_c_int + REAL(c_double) :: dnatoms + CHARACTER(LEN=100) :: error_msg + + IF (count /= 1 .AND. count /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'gather requires "count" to be 1 or 3 [Fortran/gather]') + END IF + + dnatoms = lmp_get_natoms(self) + IF (dnatoms > HUGE(1_c_int)) THEN + WRITE(error_msg,'(A,1X,I0,1X,A)') & + 'Cannot use library function gather with more than', & + HUGE(0_c_int), 'atoms [Fortran/gather]' + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, error_msg) + END IF + natoms = NINT(dnatoms, c_int) + + Cname = f2c_string(name) + IF (ALLOCATED(data)) DEALLOCATE(data) + ALLOCATE(data(natoms*count)) + Cdata = C_LOC(data(1)) + CALL lammps_gather(self%handle, Cname, Ctype, count, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_gather_double + + ! equivalent function to lammps_gather_concat (for ints) + SUBROUTINE lmp_gather_concat_int(self, name, count, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), INTENT(IN) :: count + INTEGER(c_int), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data + TYPE(c_ptr) :: Cdata, Cname + INTEGER(c_int) :: natoms + INTEGER(c_int), PARAMETER :: Ctype = 0_c_int + REAL(c_double) :: dnatoms + CHARACTER(LEN=100) :: error_msg + + IF (count /= 1 .AND. count /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'gather_concat requires "count" to be 1 or 3 [Fortran/gather_concat]') + END IF + + dnatoms = lmp_get_natoms(self) + IF (dnatoms > HUGE(1_c_int)) THEN + WRITE(error_msg,'(A,1X,I0,1X,A)') & + 'Cannot use library function gather_concat with more than', & + HUGE(0_c_int), 'atoms [Fortran/gather_concat]' + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, error_msg) + END IF + natoms = NINT(dnatoms, c_int) + + Cname = f2c_string(name) + IF (ALLOCATED(data)) DEALLOCATE(data) + ALLOCATE(data(natoms*count)) + Cdata = C_LOC(data(1)) + CALL lammps_gather_concat(self%handle, Cname, Ctype, count, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_gather_concat_int + + ! equivalent function to lammps_gather_concat (for doubles) + SUBROUTINE lmp_gather_concat_double(self, name, count, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), INTENT(IN) :: count + REAL(c_double), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data + TYPE(c_ptr) :: Cdata, Cname + INTEGER(c_int) :: natoms + INTEGER(c_int), PARAMETER :: Ctype = 1_c_int + REAL(c_double) :: dnatoms + CHARACTER(LEN=100) :: error_msg + + IF (count /= 1 .AND. count /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'gather_concat requires "count" to be 1 or 3 [Fortran/gather_concat]') + END IF + + dnatoms = lmp_get_natoms(self) + IF (dnatoms > HUGE(1_c_int)) THEN + WRITE(error_msg,'(A,1X,I0,1X,A)') & + 'Cannot use library function gather_concat with more than', & + HUGE(0_c_int), 'atoms [Fortran/gather_concat]' + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, error_msg) + END IF + natoms = NINT(dnatoms, c_int) + + Cname = f2c_string(name) + IF (ALLOCATED(data)) DEALLOCATE(data) + ALLOCATE(data(natoms*count)) + Cdata = C_LOC(data(1)) + CALL lammps_gather_concat(self%handle, Cname, Ctype, count, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_gather_concat_double + + ! equivalent function to lammps_gather_subset (for integers) + SUBROUTINE lmp_gather_subset_int(self, name, count, ids, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), INTENT(IN) :: count + INTEGER(c_int), DIMENSION(:), TARGET, INTENT(IN) :: ids + INTEGER(c_int), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data + INTEGER(c_int) :: ndata + TYPE(c_ptr) :: Cdata, Cname, Cids + INTEGER(c_int), PARAMETER :: Ctype = 0_c_int + + IF (count /= 1 .AND. count /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'gather_subset requires "count" to be 1 or 3 [Fortran/gather_subset]') + END IF + + ndata = SIZE(ids, KIND=c_int) + + Cname = f2c_string(name) + IF (ALLOCATED(data)) DEALLOCATE(data) + ALLOCATE(data(ndata*count)) + data = -1_c_int + Cdata = C_LOC(data(1)) + Cids = C_LOC(ids(1)) + CALL lammps_gather_subset(self%handle, Cname, Ctype, count, & + ndata, Cids, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_gather_subset_int + + ! equivalent function to lammps_gather_subset (for doubles) + SUBROUTINE lmp_gather_subset_double(self, name, count, ids, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), INTENT(IN) :: count + INTEGER(c_int), DIMENSION(:), TARGET, INTENT(IN) :: ids + REAL(c_double), DIMENSION(:), ALLOCATABLE, TARGET, INTENT(OUT) :: data + INTEGER(c_int) :: ndata + TYPE(c_ptr) :: Cdata, Cname, Cids + INTEGER(c_int), PARAMETER :: Ctype = 1_c_int + + IF (count /= 1 .AND. count /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'gather_subset requires "count" to be 1 or 3 [Fortran/gather_subset]') + END IF + + ndata = SIZE(ids, KIND=c_int) + + Cname = f2c_string(name) + IF (ALLOCATED(data)) DEALLOCATE(data) + ALLOCATE(data(ndata*count)) + Cdata = C_LOC(data(1)) + Cids = C_LOC(ids(1)) + CALL lammps_gather_subset(self%handle, Cname, Ctype, count, & + ndata, Cids, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_gather_subset_double + + ! equivalent function to lammps_scatter (for integers) + SUBROUTINE lmp_scatter_int(self, name, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), DIMENSION(:), TARGET :: data + INTEGER(c_int) :: natoms, Ccount + INTEGER(c_int), PARAMETER :: Ctype = 0_c_int + TYPE(c_ptr) :: Cname, Cdata + REAL(c_double) :: dnatoms + CHARACTER(LEN=100) :: error_msg + + dnatoms = lmp_get_natoms(self) + IF (dnatoms > HUGE(1_c_int)) THEN + WRITE(error_msg,'(A,1X,I0,1X,A)') & + 'Cannot use library function scatter with more than', & + HUGE(0_c_int), 'atoms [Fortran/scatter]' + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, error_msg) + END IF + natoms = NINT(dnatoms, c_int) + + Cname = f2c_string(name) + Cdata = C_LOC(data(1)) + Ccount = SIZE(data) / natoms + + IF (Ccount /= 1 .AND. Ccount /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'lammps_scatter requires either 1 or 3 data per atom') + END IF + CALL lammps_scatter(self%handle, Cname, Ctype, Ccount, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_scatter_int + + ! equivalent function to lammps_scatter (for doubles) + SUBROUTINE lmp_scatter_double(self, name, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + REAL(c_double), DIMENSION(:), TARGET :: data + INTEGER(c_int) :: natoms, Ccount + INTEGER(c_int), PARAMETER :: Ctype = 1_c_int + TYPE(c_ptr) :: Cname, Cdata + REAL(c_double) :: dnatoms + CHARACTER(LEN=100) :: error_msg + + dnatoms = lmp_get_natoms(self) + IF (dnatoms > HUGE(1_c_int)) THEN + WRITE(error_msg,'(A,1X,I0,1X,A)') & + 'Cannot use library function scatter with more than', & + HUGE(0_c_int), 'atoms [Fortran/scatter]' + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, error_msg) + END IF + natoms = NINT(dnatoms, c_int) + + Cname = f2c_string(name) + Cdata = C_LOC(data(1)) + Ccount = SIZE(data) / natoms + + IF (Ccount /= 1 .AND. Ccount /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'scatter requires either 1 or 3 data per atom [Fortran/scatter]') + END IF + CALL lammps_scatter(self%handle, Cname, Ctype, Ccount, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_scatter_double + + ! equivalent function to lammps_scatter_subset (for integers) + SUBROUTINE lmp_scatter_subset_int(self, name, ids, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), DIMENSION(:), TARGET :: ids + INTEGER(c_int), DIMENSION(:), TARGET :: data + INTEGER(c_int), PARAMETER :: Ctype = 0_c_int + INTEGER(c_int) :: Cndata, Ccount + TYPE(c_ptr) :: Cdata, Cname, Cids + + Cndata = SIZE(ids, KIND=c_int) + Ccount = SIZE(data, KIND=c_int) / Cndata + IF (Ccount /= 1 .AND. Ccount /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'scatter_subset requires either 1 or 3 data per atom & + &[Fortran/scatter_subset]') + END IF + + Cname = f2c_string(name) + Cdata = C_LOC(data(1)) + Cids = C_LOC(ids(1)) + CALL lammps_scatter_subset(self%handle, Cname, Ctype, Ccount, & + Cndata, Cids, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_scatter_subset_int + + ! equivalent function to lammps_scatter_subset (for doubles) + SUBROUTINE lmp_scatter_subset_double(self, name, ids, data) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER(c_int), DIMENSION(:), TARGET :: ids + REAL(c_double), DIMENSION(:), TARGET :: data + INTEGER(c_int), PARAMETER :: Ctype = 1_c_int + INTEGER(c_int) :: Cndata, Ccount + TYPE(c_ptr) :: Cdata, Cname, Cids + + Cndata = SIZE(ids, KIND=c_int) + Ccount = SIZE(data, KIND=c_int) / Cndata + IF (Ccount /= 1 .AND. Ccount /= 3) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'scatter_subset requires either 1 or 3 data per atom') + END IF + + Cname = f2c_string(name) + Cdata = C_LOC(data(1)) + Cids = C_LOC(ids(1)) + CALL lammps_scatter_subset(self%handle, Cname, Ctype, Ccount, & + Cndata, Cids, Cdata) + CALL lammps_free(Cname) + END SUBROUTINE lmp_scatter_subset_double + ! equivalent function to lammps_create_atoms (int ids or id absent) SUBROUTINE lmp_create_atoms_int(self, id, type, x, v, image, bexpand) CLASS(lammps), INTENT(IN) :: self @@ -1951,7 +2406,12 @@ CONTAINS CALL lammps_neighlist_element_neighbors(self%handle, idx, element, iatom, & numneigh, Cneighbors) - CALL C_F_POINTER(Cneighbors, neighbors, [numneigh]) + IF (C_ASSOCIATED(Cneighbors)) THEN + CALL C_F_POINTER(Cneighbors, neighbors, [numneigh]) + ELSE + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'Pointer returned from lammps_neighlist_element_neighbors is NULL') + END IF END SUBROUTINE lmp_neighlist_element_neighbors ! equivalent function to lammps_version @@ -2326,6 +2786,306 @@ CONTAINS END IF END SUBROUTINE lmp_decode_image_flags_bigbig + ! equivalent function to lammps_set_fix_external_callback for -DSMALLSMALL + ! note that "caller" is wrapped into a fix_external_data derived type along + ! with the fix id and the Fortran calling function. + SUBROUTINE lmp_set_fix_external_callback(self, id, callback, caller) + CLASS(lammps), TARGET, INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + EXTERNAL :: callback + CLASS(*), INTENT(IN), TARGET, OPTIONAL :: caller + TYPE(c_ptr) :: c_id, c_caller + TYPE(c_funptr) :: c_callback + INTEGER :: i, this_fix + + c_id = f2c_string(id) + IF (ALLOCATED(ext_data)) THEN + this_fix = SIZE(ext_data) + 1 + DO i = 1, SIZE(ext_data) + IF (ext_data(i)%id == id) THEN + this_fix = i + EXIT + END IF + END DO + IF (this_fix > SIZE(ext_data)) THEN + ! reallocates ext_data; this requires us to re-bind "caller" on the C + ! side to the new data structure, which likely moved to a new address + ext_data = [ext_data, fix_external_data()] ! extends ext_data by 1 + CALL rebind_external_callback_data() + END IF + ELSE + ALLOCATE(ext_data(1)) + this_fix = 1 + END IF + ext_data(this_fix)%id = id + ext_data(this_fix)%lammps_instance => self + + IF (SIZE_TAGINT == 4_c_int .AND. SIZE_BIGINT == 4_c_int) THEN + ! -DSMALLSMALL + c_callback = C_FUNLOC(callback_wrapper_smallsmall) + CALL set_fix_external_callback_smallsmall(this_fix, callback) + ELSE IF (SIZE_TAGINT == 8_c_int .AND. SIZE_BIGINT == 8_c_int) THEN + ! -DBIGBIG + c_callback = C_FUNLOC(callback_wrapper_bigbig) + CALL set_fix_external_callback_bigbig(this_fix, callback) + ELSE + ! -DSMALLBIG + c_callback = C_FUNLOC(callback_wrapper_smallbig) + CALL set_fix_external_callback_smallbig(this_fix, callback) + END IF + + IF (PRESENT(caller)) THEN + ext_data(this_fix)%caller => caller + ELSE + NULLIFY(ext_data(this_fix)%caller) + END IF + c_caller = C_LOC(ext_data(this_fix)) + CALL lammps_set_fix_external_callback(self%handle, c_id, c_callback, & + c_caller) + CALL lammps_free(c_id) + END SUBROUTINE lmp_set_fix_external_callback + + ! Wrappers to assign callback pointers with explicit interfaces + SUBROUTINE set_fix_external_callback_smallsmall(id, callback) + INTEGER, INTENT(IN) :: id + PROCEDURE(external_callback_smallsmall) :: callback + + ext_data(id)%callback_smallsmall => callback + END SUBROUTINE set_fix_external_callback_smallsmall + + SUBROUTINE set_fix_external_callback_smallbig(id, callback) + INTEGER, INTENT(IN) :: id + PROCEDURE(external_callback_smallbig) :: callback + + ext_data(id)%callback_smallbig => callback + END SUBROUTINE set_fix_external_callback_smallbig + + SUBROUTINE set_fix_external_callback_bigbig(id, callback) + INTEGER, INTENT(IN) :: id + PROCEDURE(external_callback_bigbig) :: callback + + ext_data(id)%callback_bigbig => callback + END SUBROUTINE set_fix_external_callback_bigbig + + ! subroutine that re-binds all external callback data after a reallocation + SUBROUTINE rebind_external_callback_data() + INTEGER :: i + TYPE(c_ptr) :: c_id, c_caller + TYPE(c_funptr) :: c_callback + + DO i = 1, SIZE(ext_data) - 1 + c_id = f2c_string(ext_data(i)%id) + c_caller = C_LOC(ext_data(i)) + IF (SIZE_TAGINT == 4_c_int .AND. SIZE_BIGINT == 4_c_int) THEN + c_callback = C_FUNLOC(callback_wrapper_smallsmall) + ELSE IF (SIZE_TAGINT == 8_c_int .AND. SIZE_BIGINT == 8_c_int) THEN + c_callback = C_FUNLOC(callback_wrapper_bigbig) + ELSE + c_callback = C_FUNLOC(callback_wrapper_smallbig) + END IF + CALL lammps_set_fix_external_callback( & + ext_data(i)%lammps_instance%handle, c_id, c_callback, c_caller) + CALL lammps_free(c_id) + END DO + END SUBROUTINE rebind_external_callback_data + + ! companions to lmp_set_fix_external_callback to change interface + SUBROUTINE callback_wrapper_smallsmall(caller, timestep, nlocal, ids, x, & + fexternal) BIND(C) + TYPE(c_ptr), INTENT(IN), VALUE :: caller + INTEGER(c_int), INTENT(IN), VALUE :: timestep + INTEGER(c_int), INTENT(IN), VALUE :: nlocal + TYPE(c_ptr), INTENT(IN), VALUE :: ids, x, fexternal + TYPE(c_ptr), DIMENSION(:), POINTER :: x0, f0 + INTEGER(c_int), DIMENSION(:), POINTER :: f_ids => NULL() + REAL(c_double), DIMENSION(:,:), POINTER :: f_x => NULL(), & + f_fexternal => NULL() + TYPE(fix_external_data), POINTER :: f_caller => NULL() + + CALL C_F_POINTER(ids, f_ids, [nlocal]) + CALL C_F_POINTER(x, x0, [nlocal]) + CALL C_F_POINTER(x0(1), f_x, [3, nlocal]) + CALL C_F_POINTER(fexternal, f0, [nlocal]) + CALL C_F_POINTER(f0(1), f_fexternal, [3, nlocal]) + IF (C_ASSOCIATED(caller)) THEN + CALL C_F_POINTER(caller, f_caller) + CALL f_caller%callback_smallsmall(f_caller%caller, timestep, f_ids, & + f_x, f_fexternal) + ELSE + CALL lmp_error(f_caller%lammps_instance, & + LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'Got null pointer from "caller"; this should never happen;& + & please report a bug') + END IF + END SUBROUTINE callback_wrapper_smallsmall + + SUBROUTINE callback_wrapper_smallbig(caller, timestep, nlocal, ids, x, & + fexternal) BIND(C) + TYPE(c_ptr), INTENT(IN), VALUE :: caller + INTEGER(c_int64_t), INTENT(IN), VALUE :: timestep + INTEGER(c_int), INTENT(IN), VALUE :: nlocal + TYPE(c_ptr), INTENT(IN), VALUE :: ids, x, fexternal + TYPE(c_ptr), DIMENSION(:), POINTER :: x0, f0 + INTEGER(c_int), DIMENSION(:), POINTER :: f_ids => NULL() + REAL(c_double), DIMENSION(:,:), POINTER :: f_x => NULL(), & + f_fexternal => NULL() + TYPE(fix_external_data), POINTER :: f_caller => NULL() + + CALL C_F_POINTER(ids, f_ids, [nlocal]) + CALL C_F_POINTER(x, x0, [nlocal]) + CALL C_F_POINTER(x0(1), f_x, [3, nlocal]) + CALL C_F_POINTER(fexternal, f0, [nlocal]) + CALL C_F_POINTER(f0(1), f_fexternal, [3, nlocal]) + IF (C_ASSOCIATED(caller)) THEN + CALL C_F_POINTER(caller, f_caller) + CALL f_caller%callback_smallbig(f_caller%caller, timestep, f_ids, f_x, & + f_fexternal) + ELSE + CALL lmp_error(f_caller%lammps_instance, & + LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'Got null pointer from "caller"; this should never happen;& + & please report a bug') + END IF + END SUBROUTINE callback_wrapper_smallbig + + SUBROUTINE callback_wrapper_bigbig(caller, timestep, nlocal, ids, x, & + fexternal) BIND(C) + TYPE(c_ptr), INTENT(IN), VALUE :: caller + INTEGER(c_int64_t), INTENT(IN), VALUE :: timestep + INTEGER(c_int), INTENT(IN), VALUE :: nlocal + TYPE(c_ptr), INTENT(IN), VALUE :: ids, x, fexternal + TYPE(c_ptr), DIMENSION(:), POINTER :: x0, f0 + INTEGER(c_int64_t), DIMENSION(:), POINTER :: f_ids => NULL() + REAL(c_double), DIMENSION(:,:), POINTER :: f_x => NULL(), & + f_fexternal => NULL() + TYPE(fix_external_data), POINTER :: f_caller => NULL() + + CALL C_F_POINTER(ids, f_ids, [nlocal]) + CALL C_F_POINTER(x, x0, [nlocal]) + CALL C_F_POINTER(x0(1), f_x, [3, nlocal]) + CALL C_F_POINTER(fexternal, f0, [nlocal]) + CALL C_F_POINTER(f0(1), f_fexternal, [3, nlocal]) + IF (C_ASSOCIATED(caller)) THEN + CALL C_F_POINTER(caller, f_caller) + CALL f_caller%callback_bigbig(f_caller%caller, timestep, f_ids, f_x, & + f_fexternal) + ELSE + CALL lmp_error(f_caller%lammps_instance, & + LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'Got null pointer from "caller"; this should never happen;& + & please report a bug') + END IF + END SUBROUTINE callback_wrapper_bigbig + + ! equivalent function to lammps_fix_external_get_force + FUNCTION lmp_fix_external_get_force(self, id) RESULT(fexternal) + CLASS(lammps), TARGET, INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + TYPE(lammps_fix_data) :: fexternal + TYPE(c_ptr) :: ptr, Cid + TYPE(c_ptr), DIMENSION(:), POINTER :: f + INTEGER(c_int) :: nmax + + Cid = f2c_string(id) + ptr = lammps_fix_external_get_force(self%handle, Cid) + nmax = lmp_extract_setting(self, 'nmax') + CALL C_F_POINTER(ptr, f, [nmax]) + fexternal%datatype = DATA_DOUBLE_2D + fexternal%lammps_instance => self + CALL C_F_POINTER(f(1), fexternal%r64_mat, [3, nmax]) + CALL lammps_free(Cid) + END FUNCTION lmp_fix_external_get_force + + SUBROUTINE lmp_fix_external_set_energy_global(self, id, eng) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + REAL(c_double), INTENT(IN) :: eng + TYPE(c_ptr) :: Cid + + Cid = f2c_string(id) + CALL lammps_fix_external_set_energy_global(self%handle, Cid, eng) + CALL lammps_free(Cid) + END SUBROUTINE lmp_fix_external_set_energy_global + + SUBROUTINE lmp_fix_external_set_virial_global(self, id, virial) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + REAL(c_double), DIMENSION(6), TARGET, INTENT(IN) :: virial + TYPE(c_ptr) :: Cid, Cvirial + + Cid = f2c_string(id) + Cvirial = C_LOC(virial(1)) + CALL lammps_fix_external_set_virial_global(self%handle, Cid, Cvirial) + CALL lammps_free(Cid) + END SUBROUTINE lmp_fix_external_set_virial_global + + SUBROUTINE lmp_fix_external_set_energy_peratom(self, id, eng) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + REAL(c_double), DIMENSION(:), TARGET, INTENT(IN) :: eng + TYPE(c_ptr) :: Cid, Ceng + INTEGER(c_int) :: nlocal + + nlocal = lmp_extract_setting(self, 'nlocal') + IF (SIZE(eng) < nlocal) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'Array "eng" should be length nlocal or greater & + &[Fortran/fix_external_set_energy_peratom]') + END IF + Cid = f2c_string(id) + Ceng = C_LOC(eng) + CALL lammps_fix_external_set_energy_peratom(self%handle, Cid, Ceng) + CALL lammps_free(Cid) + END SUBROUTINE lmp_fix_external_set_energy_peratom + + SUBROUTINE lmp_fix_external_set_virial_peratom(self, id, virial) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + REAL(c_double), DIMENSION(:,:), TARGET, INTENT(IN) :: virial + TYPE(c_ptr) :: Cid, Cvirial + TYPE(c_ptr), DIMENSION(:), ALLOCATABLE, TARGET :: Cptr + INTEGER(c_int) :: nlocal, i + + nlocal = lmp_extract_setting(self, 'nlocal') + IF (SIZE(virial,2) < nlocal .OR. SIZE(virial,1) /= 6) THEN + CALL lmp_error(self, LMP_ERROR_ALL + LMP_ERROR_WORLD, & + 'Array "virial" should be size 6 x nlocal or greater & + &[Fortran/fix_external_set_energy_peratom]') + END IF + Cid = f2c_string(id) + ALLOCATE(Cptr(nlocal)) + DO i = 1, nlocal + Cptr(i) = C_LOC(virial(1,i)) + END DO + Cvirial = C_LOC(Cptr(1)) + CALL lammps_fix_external_set_virial_peratom(self%handle, Cid, Cvirial) + CALL lammps_free(Cid) + DEALLOCATE(Cptr) + END SUBROUTINE lmp_fix_external_set_virial_peratom + + SUBROUTINE lmp_fix_external_set_vector_length(self, id, length) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + INTEGER(c_int), INTENT(IN) :: length + TYPE(c_ptr) :: Cid + + Cid = f2c_string(id) + CALL lammps_fix_external_set_vector_length(self%handle, Cid, length) + CALL lammps_free(Cid) + END SUBROUTINE lmp_fix_external_set_vector_length + + SUBROUTINE lmp_fix_external_set_vector(self, id, idx, val) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: id + INTEGER(c_int), INTENT(IN) :: idx + REAL(c_double), INTENT(IN) :: val + TYPE(c_ptr) :: Cid + + Cid = f2c_string(id) + CALL lammps_fix_external_set_vector(self%handle, Cid, idx, val) + CALL lammps_free(Cid) + END SUBROUTINE lmp_fix_external_set_vector + ! equivalent function to lammps_flush_buffers SUBROUTINE lmp_flush_buffers(self) CLASS(lammps), INTENT(IN) :: self diff --git a/src/library.cpp b/src/library.cpp index 1e3c96cdcf..2c16cc6615 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -685,7 +685,7 @@ void lammps_commands_string(void *handle, const char *str) \verbatim embed:rst This number may be very large when running large simulations across -multiple processors. Depending on compile time choices, LAMMPS may be +multiple processes. Depending on compile time choices, LAMMPS may be using either 32-bit or a 64-bit integer to store this number. For portability this function returns thus a double precision floating point number, which can represent up to a 53-bit signed @@ -2263,7 +2263,7 @@ int lammps_set_variable(void *handle, char *name, char *str) // Library functions for scatter/gather operations of data // ---------------------------------------------------------------------- -/** Gather the named atom-based entity for all atoms across all processors, +/** Gather the named atom-based entity for all atoms across all processes, * in order. * \verbatim embed:rst @@ -2282,6 +2282,8 @@ x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], :math:`\dots`); *natoms*), as queried by :cpp:func:`lammps_get_natoms`, :cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`. +This function is not compatible with ``-DLAMMPS_BIGBIG``. + \endverbatim * * \param handle pointer to a previously created LAMMPS instance @@ -2301,7 +2303,8 @@ x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], :math:`\dots`); Allreduce to sum vector into data across all procs ------------------------------------------------------------------------- */ -void lammps_gather_atoms(void *handle, char *name, int type, int count, void *data) +void lammps_gather_atoms(void *handle, const char *name, int type, int count, + void *data) { auto lmp = (LAMMPS *) handle; @@ -2415,7 +2418,7 @@ void lammps_gather_atoms(void *handle, char *name, int type, int count, void *da END_CAPTURE } -/** Gather the named atom-based entity for all atoms across all processors, +/** Gather the named atom-based entity for all atoms across all processes, * unordered. * \verbatim embed:rst @@ -2432,12 +2435,14 @@ of atoms, use :cpp:func:`lammps_gather_atoms_subset`. The *data* array will be in groups of *count* values, with *natoms* groups total, but not in order by atom ID (e.g., if *name* is *x* and *count* -is 3, then *data* might be something like = x[10][0], x[10][1], x[10][2], +is 3, then *data* might be something like x[10][0], x[10][1], x[10][2], x[2][0], x[2][1], x[2][2], x[4][0], :math:`\dots`); *data* must be pre-allocated by the caller to length (*count* :math:`\times` *natoms*), as queried by :cpp:func:`lammps_get_natoms`, :cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`. +This function is not compatible with ``-DLAMMPS_BIGBIG``. + \endverbatim * * \param handle: pointer to a previously created LAMMPS instance @@ -2456,7 +2461,8 @@ queried by :cpp:func:`lammps_get_natoms`, Allgather Nlocal atoms from each proc into data ------------------------------------------------------------------------- */ -void lammps_gather_atoms_concat(void *handle, char *name, int type, int count, void *data) +void lammps_gather_atoms_concat(void *handle, const char *name, int type, + int count, void *data) { auto lmp = (LAMMPS *) handle; @@ -2599,6 +2605,8 @@ x[100][2], x[57][0], x[57][1], x[57][2], x[210][0], :math:`\dots`); *data* must be pre-allocated by the caller to length (*count* :math:`\times` *ndata*). +This function is not compatible with ``-DLAMMPS_BIGBIG``. + \endverbatim * * \param handle: pointer to a previously created LAMMPS instance @@ -2621,8 +2629,8 @@ x[100][2], x[57][0], x[57][1], x[57][2], x[210][0], :math:`\dots`); Allreduce to sum vector into data across all procs ------------------------------------------------------------------------- */ -void lammps_gather_atoms_subset(void *handle, char *name, int type, int count, - int ndata, int *ids, void *data) +void lammps_gather_atoms_subset(void *handle, const char *name, int type, + int count, int ndata, int *ids, void *data) { auto lmp = (LAMMPS *) handle; @@ -2745,20 +2753,22 @@ void lammps_gather_atoms_subset(void *handle, char *name, int type, int count, END_CAPTURE } -/** Scatter the named atom-based entities in *data* to all processors. +/** Scatter the named atom-based entities in *data* to all processes. * \verbatim embed:rst This subroutine takes data stored in a one-dimensional array supplied by the -user and scatters them to all atoms on all processors. The data must be +user and scatters them to all atoms on all processes. The data must be ordered by atom ID, with the requirement that the IDs be consecutive. Use :cpp:func:`lammps_scatter_atoms_subset` to scatter data for some (or all) atoms, unordered. The *data* array needs to be ordered in groups of *count* values, sorted by atom ID (e.g., if *name* is *x* and *count* = 3, then -*data* = x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], -:math:`\dots`); *data* must be of length (*count* :math:`\times` *natoms*). +*data* = {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], +:math:`\dots`}); *data* must be of length (*count* :math:`\times` *natoms*). + +This function is not compatible with ``-DLAMMPS_BIGBIG``. \endverbatim * @@ -2768,7 +2778,7 @@ atom ID (e.g., if *name* is *x* and *count* = 3, then * \param count number of per-atom values (e.g., 1 for *type* or *charge*, * 3 for *x* or *f*); use *count* = 3 with *image* if you have * a single image flag packed into (*x*,*y*,*z*) components. - * \param data per-atom values packed in a 1-dimensional array of length + * \param data per-atom values packed in a one-dimensional array of length * *natoms* \* *count*. * */ @@ -2778,7 +2788,8 @@ atom ID (e.g., if *name* is *x* and *count* = 3, then loop over Natoms, if I own atom ID, set its values from data ------------------------------------------------------------------------- */ -void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *data) +void lammps_scatter_atoms(void *handle, const char *name, int type, int count, + void *data) { auto lmp = (LAMMPS *) handle; @@ -2879,12 +2890,12 @@ void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *d } /** Scatter the named atom-based entities in *data* from a subset of atoms - * to all processors. + * to all processes. * \verbatim embed:rst This subroutine takes data stored in a one-dimensional array supplied by the -user and scatters them to a subset of atoms on all processors. The array +user and scatters them to a subset of atoms on all processes. The array *data* contains data associated with atom IDs, but there is no requirement that the IDs be consecutive, as they are provided in a separate array. Use :cpp:func:`lammps_scatter_atoms` to scatter data for all atoms, in order. @@ -2895,6 +2906,8 @@ to be the array {x[1][0], x[1][1], x[1][2], x[100][0], x[100][1], x[100][2], x[57][0], x[57][1], x[57][2]}, then *count* = 3, *ndata* = 3, and *ids* would be {1, 100, 57}. +This function is not compatible with ``-DLAMMPS_BIGBIG``. + \endverbatim * * \param handle: pointer to a previously created LAMMPS instance @@ -2928,8 +2941,8 @@ be {1, 100, 57}. loop over Ndata, if I own atom ID, set its values from data ------------------------------------------------------------------------- */ -void lammps_scatter_atoms_subset(void *handle, char *name, int type, int count, - int ndata, int *ids, void *data) +void lammps_scatter_atoms_subset(void *handle, const char *name, int type, + int count, int ndata, int *ids, void *data) { auto lmp = (LAMMPS *) handle; @@ -3146,6 +3159,43 @@ void lammps_gather_bonds(void *handle, void *data) END_CAPTURE } +/** Gather the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based entities + * from all processes, in order by atom ID. + * +\verbatim embed:rst + +This subroutine gathers data from all processes and stores them in a one-dimensional array +allocated by the user. The array *data* will be ordered by atom ID, which requires consecutive IDs +(1 to *natoms*\ ). If you need a similar array but for non-consecutive atom IDs, see +:cpp:func:`lammps_gather_concat`; for a similar array but for a subset of atoms, see +:cpp:func:`lammps_gather_subset`. + +The *data* array will be ordered in groups of *count* values, sorted by atom ID (e.g., if *name* is +*x*, then *data* is {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], +:math:`\dots`}); *data* must be pre-allocated by the caller to the correct length +(*count*\ :math:`{}\times{}`\ *natoms*), as queried by :cpp:func:`lammps_get_natoms`, +:cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`. + +This function will return an error if fix or compute data are requested and the fix or compute ID +given does not have per-atom data. + +This function is not compatible with ``-DLAMMPS_BIGBIG``. + +\endverbatim + * + * \param handle pointer to a previously created LAMMPS instance + * \param name desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix + * data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix + * property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix + * property/atom vectors with *count* > 1) + * \param type 0 for ``int`` values, 1 for ``double`` values + * \param count number of per-atom values (e.g., 1 for *type* or *charge*, 3 for *x* or *f*); + * use *count* = 3 with *image* if you want the image flags unpacked into + * (*x*,*y*,*z*) components. + * \param data per-atom values packed into a one-dimensional array of length + * *natoms* \* *count*. + * + */ /* ---------------------------------------------------------------------- Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France) gather the named atom-based entity for all atoms @@ -3172,7 +3222,7 @@ void lammps_gather_bonds(void *handle, void *data) Allreduce to sum vector into data across all procs ------------------------------------------------------------------------- */ -void lammps_gather(void *handle, char *name, int type, int count, void *data) +void lammps_gather(void *handle, const char *name, int type, int count, void *data) { auto lmp = (LAMMPS *) handle; @@ -3381,6 +3431,44 @@ void lammps_gather(void *handle, char *name, int type, int count, void *data) END_CAPTURE } +/** Gather the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based entities + * from all processes, unordered. + * +\verbatim embed:rst + +This subroutine gathers data for all atoms and stores them in a one-dimensional array allocated by +the user. The data will be a concatenation of chunks from each processor's owned atoms, in +whatever order the atoms are in on each processor. This process has no requirement that the atom +IDs be consecutive. If you need the ID of each atom, you can do another call to either +:cpp:func:`lammps_gather_atoms_concat` or :cpp:func:`lammps_gather_concat` with *name* set to +``id``. If you have consecutive IDs and want the data to be in order, use +:cpp:func:`lammps_gather`; for a similar array but for a subset of atoms, use +:cpp:func:`lammps_gather_subset`. + +The *data* array will be in groups of *count* values, with *natoms* groups total, but not in order +by atom ID (e.g., if *name* is *x* and *count* is 3, then *data* might be something like +{x[10][0], x[10][1], x[10][2], x[2][0], x[2][1], x[2][2], x[4][0], :math:`\dots`}); *data* must be +pre-allocated by the caller to length (*count* :math:`\times` *natoms*), as queried by +:cpp:func:`lammps_get_natoms`, :cpp:func:`lammps_extract_global`, or +:cpp:func:`lammps_extract_setting`. + +This function is not compatible with ``-DLAMMPS_BIGBIG``. + +\endverbatim + * + * \param handle: pointer to a previously created LAMMPS instance + * \param name: desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix + * data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix + * property/atom vectors with count = 1, "d2_name" or "i2_name" for fix + * property/atom vectors with count > 1) + * \param type: 0 for ``int`` values, 1 for ``double`` values + * \param count: number of per-atom values (e.g., 1 for *type* or *charge*, 3 for *x* or *f*); + * use *count* = 3 with *image* if you want the image flags unpacked into + * (*x*,*y*,*z*) components. + * \param data: per-atom values packed into a one-dimensional array of length + * *natoms* \* *count*. + * + */ /* ---------------------------------------------------------------------- Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France) gather the named atom-based entity for all atoms @@ -3407,7 +3495,8 @@ void lammps_gather(void *handle, char *name, int type, int count, void *data) Allreduce to sum vector into data across all procs ------------------------------------------------------------------------- */ -void lammps_gather_concat(void *handle, char *name, int type, int count, void *data) +void lammps_gather_concat(void *handle, const char *name, int type, int count, + void *data) { auto lmp = (LAMMPS *) handle; @@ -3633,6 +3722,41 @@ void lammps_gather_concat(void *handle, char *name, int type, int count, void *d END_CAPTURE } +/** Gather the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based entities + * from all processes for a subset of atoms. + * +\verbatim embed:rst + +This subroutine gathers data for the requested atom IDs and stores them in a one-dimensional array +allocated by the user. The data will be ordered by atom ID, but there is no requirement that the +IDs be consecutive. If you wish to return a similar array for *all* the atoms, use +:cpp:func:`lammps_gather` or :cpp:func:`lammps_gather_concat`. + +The *data* array will be in groups of *count* values, sorted by atom ID in the same order as the +array *ids* (e.g., if *name* is *x*, *count* = 3, and *ids* is {100, 57, 210}, then *data* might +look like {x[100][0], x[100][1], x[100][2], x[57][0], x[57][1], x[57][2], x[210][0], +:math:`\dots`}); *ids* must be provided by the user with length *ndata*, and *data* must be +pre-allocated by the caller to length (*count*\ :math:`{}\times{}`\ *ndata*). + +This function is not compatible with ``-DLAMMPS_BIGBIG``. + +\endverbatim + * + * \param handle: pointer to a previously created LAMMPS instance + * \param name desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix + * data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix + * property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix + * property/atom vectors with *count* > 1) + * \param type 0 for ``int`` values, 1 for ``double`` values + * \param count number of per-atom values (e.g., 1 for *type* or *charge*, 3 for *x* or *f*); + * use *count* = 3 with *image* if you want the image flags unpacked into + * (*x*,*y*,*z*) components. + * \param ndata: number of atoms for which to return data (can be all of them) + * \param ids: list of *ndata* atom IDs for which to return data + * \param data per-atom values packed into a one-dimensional array of length + * *ndata* \* *count*. + * + */ /* ---------------------------------------------------------------------- Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France) gather the named atom-based entity for all atoms @@ -3659,9 +3783,8 @@ void lammps_gather_concat(void *handle, char *name, int type, int count, void *d Allreduce to sum vector into data across all procs ------------------------------------------------------------------------- */ -void lammps_gather_subset(void *handle, char *name, - int type, int count, - int ndata, int *ids, void *data) +void lammps_gather_subset(void *handle, const char *name, int type, int count, + int ndata, int *ids, void *data) { auto lmp = (LAMMPS *) handle; @@ -3884,6 +4007,37 @@ void lammps_gather_subset(void *handle, char *name, END_CAPTURE } +/** Scatter the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based + * entity in *data* to all processes. + * +\verbatim embed:rst + +This subroutine takes data stored in a one-dimensional array supplied by the user and scatters +them to all atoms on all processes. The data must be ordered by atom ID, with the requirement that +the IDs be consecutive. Use :cpp:func:`lammps_scatter_subset` to scatter data for some (or all) +atoms, unordered. + +The *data* array needs to be ordered in groups of *count* values, sorted by atom ID (e.g., if +*name* is *x* and *count* = 3, then *data* = {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], +x[1][2], x[2][0], :math:`\dots`}); *data* must be of length (*count* :math:`\times` *natoms*). + +This function is not compatible with ``-DLAMMPS_BIGBIG``. + +\endverbatim + * + * \param handle pointer to a previously created LAMMPS instance + * \param name desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix + * data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix + * property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix + * property/atom vectors with *count* > 1) + * \param type 0 for ``int`` values, 1 for ``double`` values + * \param count number of per-atom values (e.g., 1 for *type* or *charge*, + * 3 for *x* or *f*); use *count* = 3 with *image* if you have + * a single image flag packed into (*x*,*y*,*z*) components. + * \param data per-atom values packed in a one-dimensional array of length + * *natoms* \* *count*. + * + */ /* ---------------------------------------------------------------------- Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France) scatter the named atom-based entity in data to all atoms @@ -3908,7 +4062,8 @@ void lammps_gather_subset(void *handle, char *name, Allreduce to sum vector into data across all procs ------------------------------------------------------------------------- */ -void lammps_scatter(void *handle, char *name, int type, int count, void *data) +void lammps_scatter(void *handle, const char *name, int type, int count, + void *data) { auto lmp = (LAMMPS *) handle; @@ -4103,6 +4258,42 @@ void lammps_scatter(void *handle, char *name, int type, int count, void *data) END_CAPTURE } +/** Scatter the named per-atom, per-atom fix, per-atom compute, or fix property/atom-based + * entities in *data* from a subset of atoms to all processes. + * +\verbatim embed:rst + +This subroutine takes data stored in a one-dimensional array supplied by the +user and scatters them to a subset of atoms on all processes. The array +*data* contains data associated with atom IDs, but there is no requirement that +the IDs be consecutive, as they are provided in a separate array. +Use :cpp:func:`lammps_scatter` to scatter data for all atoms, in order. + +The *data* array needs to be organized in groups of *count* values, with the +groups in the same order as the array *ids*. For example, if you want *data* +to be the array {x[1][0], x[1][1], x[1][2], x[100][0], x[100][1], x[100][2], +x[57][0], x[57][1], x[57][2]}, then *count* = 3, *ndata* = 3, and *ids* would +be {1, 100, 57}. + +This function is not compatible with ``-DLAMMPS_BIGBIG``. + +\endverbatim + * + * \param handle: pointer to a previously created LAMMPS instance + * \param name desired quantity (e.g., "x" or "f" for atom properties, "f_id" for per-atom fix + * data, "c_id" for per-atom compute data, "d_name" or "i_name" for fix + * property/atom vectors with *count* = 1, "d2_name" or "i2_name" for fix + * property/atom vectors with *count* > 1) + * \param type: 0 for ``int`` values, 1 for ``double`` values + * \param count: number of per-atom values (e.g., 1 for *type* or *charge*, + * 3 for *x* or *f*); use *count* = 3 with "image" if you want + * single image flags unpacked into (*x*,*y*,*z*) + * \param ndata: number of atoms listed in *ids* and *data* arrays + * \param ids: list of *ndata* atom IDs to scatter data to + * \param data per-atom values packed in a 1-dimensional array of length + * *ndata* \* *count*. + * + */ /* ---------------------------------------------------------------------- Contributing author: Thomas Swinburne (CNRS & CINaM, Marseille, France) scatter the named atom-based entity in data to a subset of atoms @@ -4125,7 +4316,7 @@ void lammps_scatter(void *handle, char *name, int type, int count, void *data) loop over Ndata, if I own atom ID, set its values from data ------------------------------------------------------------------------- */ -void lammps_scatter_subset(void *handle, char *name,int type, int count, +void lammps_scatter_subset(void *handle, const char *name,int type, int count, int ndata, int *ids, void *data) { auto lmp = (LAMMPS *) handle; diff --git a/src/library.h b/src/library.h index f0f8448b79..9cd2f34d73 100644 --- a/src/library.h +++ b/src/library.h @@ -181,23 +181,23 @@ int lammps_set_variable(void *, char *, char *); * Library functions for scatter/gather operations of data * ---------------------------------------------------------------------- */ -void lammps_gather_atoms(void *handle, char *name, int type, int count, void *data); -void lammps_gather_atoms_concat(void *handle, char *name, int type, int count, void *data); -void lammps_gather_atoms_subset(void *handle, char *name, int type, int count, int ndata, int *ids, - void *data); -void lammps_scatter_atoms(void *handle, char *name, int type, int count, void *data); -void lammps_scatter_atoms_subset(void *handle, char *name, int type, int count, int ndata, int *ids, - void *data); +void lammps_gather_atoms(void *handle, const char *name, int type, int count, void *data); +void lammps_gather_atoms_concat(void *handle, const char *name, int type, int count, void *data); +void lammps_gather_atoms_subset(void *handle, const char *name, int type, int count, int ndata, + int *ids, void *data); +void lammps_scatter_atoms(void *handle, const char *name, int type, int count, void *data); +void lammps_scatter_atoms_subset(void *handle, const char *name, int type, int count, int ndata, + int *ids, void *data); void lammps_gather_bonds(void *handle, void *data); -void lammps_gather(void *handle, char *name, int type, int count, void *data); -void lammps_gather_concat(void *handle, char *name, int type, int count, void *data); -void lammps_gather_subset(void *handle, char *name, int type, int count, int ndata, int *ids, - void *data); -void lammps_scatter(void *handle, char *name, int type, int count, void *data); -void lammps_scatter_subset(void *handle, char *name, int type, int count, int ndata, int *ids, - void *data); +void lammps_gather(void *handle, const char *name, int type, int count, void *data); +void lammps_gather_concat(void *handle, const char *name, int type, int count, void *data); +void lammps_gather_subset(void *handle, const char *name, int type, int count, int ndata, + int *ids, void *data); +void lammps_scatter(void *handle, const char *name, int type, int count, void *data); +void lammps_scatter_subset(void *handle, const char *name, int type, int count, int ndata, + int *ids, void *data); #if !defined(LAMMPS_BIGBIG) int lammps_create_atoms(void *handle, int n, const int *id, const int *type, const double *x, diff --git a/unittest/fortran/CMakeLists.txt b/unittest/fortran/CMakeLists.txt index b534835a5e..93573eb297 100644 --- a/unittest/fortran/CMakeLists.txt +++ b/unittest/fortran/CMakeLists.txt @@ -86,6 +86,10 @@ if(CMAKE_Fortran_COMPILER) target_link_libraries(test_fortran_neighlist PRIVATE flammps lammps GTest::GMockMain) add_test(NAME FortranNeighlist COMMAND test_fortran_neighlist) + add_executable(test_fortran_fixexternal wrap_fixexternal.cpp test_fortran_fixexternal.f90) + target_link_libraries(test_fortran_fixexternal PRIVATE flammps lammps GTest::GMockMain) + add_test(NAME FortranFixExternal COMMAND test_fortran_fixexternal) + else() message(STATUS "Skipping Tests for the LAMMPS Fortran Module: no Fortran compiler") endif() diff --git a/unittest/fortran/test_fortran_fixexternal.f90 b/unittest/fortran/test_fortran_fixexternal.f90 new file mode 100644 index 0000000000..d46e1ec673 --- /dev/null +++ b/unittest/fortran/test_fortran_fixexternal.f90 @@ -0,0 +1,424 @@ +MODULE ext_stuff + USE, INTRINSIC :: ISO_Fortran_ENV, ONLY : error_unit + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int, c_int64_t, c_loc + USE LIBLAMMPS + IMPLICIT NONE + + INTEGER, PARAMETER :: vec_length = 8 + REAL(c_double), SAVE :: direction = 1.0_c_double + REAL(c_double), DIMENSION(:,:), POINTER, SAVE :: f3 => NULL(), f4 => NULL() + +CONTAINS + + SUBROUTINE f_lammps_reverse_direction() BIND(C) + direction = -direction + END SUBROUTINE f_lammps_reverse_direction + + SUBROUTINE f_callback_ss(instance, timestep, id, x, f) + CLASS(*), INTENT(INOUT) :: instance + INTEGER(c_int) :: timestep + INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f + REAL(c_double), DIMENSION(SIZE(id)) :: e + REAL(c_double), DIMENSION(6,SIZE(id)) :: v + + WHERE (id == 1) + f(1,:) = 1.0_c_double + f(2,:) = -1.0_c_double + f(3,:) = 1.25_c_double + e = 1.0_c_double + v(1,:) = 1.0_c_double + v(2,:) = 2.0_c_double + v(3,:) = -1.0_c_double + v(4,:) = -2.0_c_double + v(5,:) = 3.0_c_double + v(6,:) = -3.0_c_double + ELSEWHERE + f(1,:) = -1.0_c_double + f(2,:) = +1.0_c_double + f(3,:) = -1.25_c_double + e = 10.0_c_double + v(1,:) = 10.0_c_double + v(2,:) = 20.0_c_double + v(3,:) = -10.0_c_double + v(4,:) = -20.0_c_double + v(5,:) = 30.0_c_double + v(6,:) = -30.0_c_double + END WHERE + SELECT TYPE (instance) + CLASS IS (lammps) + CALL instance%fix_external_set_energy_peratom('ext1', e) + CALL instance%fix_external_set_virial_peratom('ext1', v) + CLASS DEFAULT + WRITE(error_unit,*) 'UMM...this should never happen.' + STOP 1 + END SELECT + END SUBROUTINE f_callback_ss + + SUBROUTINE f_callback_sb(instance, timestep, id, x, f) + CLASS(*), INTENT(INOUT) :: instance + INTEGER(c_int64_t) :: timestep + INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f + REAL(c_double), DIMENSION(SIZE(id)) :: e + REAL(c_double), DIMENSION(6,SIZE(id)) :: v + + WHERE (id == 1_c_int) + f(1,:) = 1.0_c_double + f(2,:) = -1.0_c_double + f(3,:) = 1.25_c_double + e = 1.0_c_double + v(1,:) = 1.0_c_double + v(2,:) = 2.0_c_double + v(3,:) = -1.0_c_double + v(4,:) = -2.0_c_double + v(5,:) = 3.0_c_double + v(6,:) = -3.0_c_double + ELSEWHERE + f(1,:) = -1.0_c_double + f(2,:) = +1.0_c_double + f(3,:) = -1.25_c_double + e = 10.0_c_double + v(1,:) = 10.0_c_double + v(2,:) = 20.0_c_double + v(3,:) = -10.0_c_double + v(4,:) = -20.0_c_double + v(5,:) = 30.0_c_double + v(6,:) = -30.0_c_double + END WHERE + SELECT TYPE (instance) + CLASS IS (lammps) + CALL instance%fix_external_set_energy_peratom('ext1', e) + CALL instance%fix_external_set_virial_peratom('ext1', v) + CLASS DEFAULT + WRITE(error_unit,*) 'UMM...this should never happen.' + STOP 1 + END SELECT + END SUBROUTINE f_callback_sb + + SUBROUTINE f_callback_bb(instance, timestep, id, x, f) + CLASS(*), INTENT(INOUT) :: instance + INTEGER(c_int64_t) :: timestep + INTEGER(c_int64_t), DIMENSION(:), INTENT(IN) :: id + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f + REAL(c_double), DIMENSION(SIZE(id)) :: e + REAL(c_double), DIMENSION(6,SIZE(id)) :: v + + WHERE (id == 1_c_int64_t) + f(1,:) = 1.0_c_double + f(2,:) = -1.0_c_double + f(3,:) = 1.25_c_double + e = 1.0_c_double + v(1,:) = 1.0_c_double + v(2,:) = 2.0_c_double + v(3,:) = -1.0_c_double + v(4,:) = -2.0_c_double + v(5,:) = 3.0_c_double + v(6,:) = -3.0_c_double + ELSEWHERE + f(1,:) = -1.0_c_double + f(2,:) = +1.0_c_double + f(3,:) = -1.25_c_double + e = 10.0_c_double + v(1,:) = 10.0_c_double + v(2,:) = 20.0_c_double + v(3,:) = -10.0_c_double + v(4,:) = -20.0_c_double + v(5,:) = 30.0_c_double + v(6,:) = -30.0_c_double + END WHERE + SELECT TYPE (instance) + CLASS IS (lammps) + CALL instance%fix_external_set_energy_peratom('ext1', e) + CALL instance%fix_external_set_virial_peratom('ext1', v) + CLASS DEFAULT + WRITE(error_unit,*) 'UMM...this should never happen.' + STOP 1 + END SELECT + END SUBROUTINE f_callback_bb + + SUBROUTINE f_callback2_ss(entity, timestep, id, x, f) + CLASS(*), INTENT(INOUT), target :: entity + INTEGER(c_int) :: timestep + INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f + + SELECT TYPE (entity) + TYPE IS (REAL(c_double)) + WHERE (id == 1_c_int) + f(1,:) = SIGN(1.0_c_double, entity) * 2.0_c_double + f(2,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double) + f(3,:) = SIGN(1.0_c_double, entity) * 2.5_c_double + ELSEWHERE + f(1,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double) + f(2,:) = SIGN(1.0_c_double, entity) * 2.0_c_double + f(3,:) = SIGN(1.0_c_double, entity) * (-2.5_c_double) + END WHERE + CLASS DEFAULT + WRITE(error_unit,'(A)') 'ERROR: Failed to resolve "entity" in& + & f_callback2_ss' + STOP 1 + END SELECT + END SUBROUTINE f_callback2_ss + + SUBROUTINE f_callback2_sb(entity, timestep, id, x, f) + CLASS(*), INTENT(INOUT), target :: entity + INTEGER(c_int64_t) :: timestep + INTEGER(c_int), DIMENSION(:), INTENT(IN) :: id + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f + + SELECT TYPE (entity) + TYPE IS (REAL(c_double)) + WHERE (id == 1_c_int) + f(1,:) = SIGN(1.0_c_double, entity) * 2.0_c_double + f(2,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double) + f(3,:) = SIGN(1.0_c_double, entity) * 2.5_c_double + ELSEWHERE + f(1,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double) + f(2,:) = SIGN(1.0_c_double, entity) * 2.0_c_double + f(3,:) = SIGN(1.0_c_double, entity) * (-2.5_c_double) + END WHERE + CLASS DEFAULT + WRITE(error_unit,'(A)') 'ERROR: Failed to resolve "entity" in& + & f_callback2_sb' + STOP 1 + END SELECT + END SUBROUTINE f_callback2_sb + + SUBROUTINE f_callback2_bb(entity, timestep, id, x, f) + CLASS(*), INTENT(INOUT), target :: entity + INTEGER(c_int64_t) :: timestep + INTEGER(c_int64_t), DIMENSION(:), INTENT(IN) :: id + REAL(c_double), DIMENSION(:,:), INTENT(IN) :: x + REAL(c_double), DIMENSION(:,:), INTENT(OUT) :: f + + SELECT TYPE (entity) + TYPE IS (REAL(c_double)) + WHERE (id == 1_c_int64_t) + f(1,:) = SIGN(1.0_c_double, entity) * 2.0_c_double + f(2,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double) + f(3,:) = SIGN(1.0_c_double, entity) * 2.5_c_double + ELSEWHERE + f(1,:) = SIGN(1.0_c_double, entity) * (-2.0_c_double) + f(2,:) = SIGN(1.0_c_double, entity) * 2.0_c_double + f(3,:) = SIGN(1.0_c_double, entity) * (-2.5_c_double) + END WHERE + CLASS DEFAULT + WRITE(error_unit,'(A)') 'ERROR: Failed to resolve "entity" in& + & f_callback2_sb' + STOP 1 + END SELECT + END SUBROUTINE f_callback2_bb +END MODULE ext_stuff + +FUNCTION f_lammps_with_args() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_ptr + USE liblammps + USE keepstuff, 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, INTRINSIC :: ISO_C_BINDING, ONLY: c_null_ptr + USE liblammps + USE keepstuff, ONLY: lmp + IMPLICIT NONE + + CALL lmp%close() + lmp%handle = c_null_ptr +END SUBROUTINE f_lammps_close + +SUBROUTINE f_lammps_setup_fix_external_callback() BIND(C) + USE LIBLAMMPS + USE keepstuff, ONLY : lmp, demo_input, cont_input, pair_input + USE ext_stuff, ONLY : vec_length + IMPLICIT NONE + + CALL lmp%commands_list(demo_input) + CALL lmp%commands_list(cont_input) + CALL lmp%commands_list(pair_input) + CALL lmp%command('neigh_modify exclude group all all') + CALL lmp%command('fix ext1 all external pf/callback 1 1') + CALL lmp%command('fix ext2 all external pf/callback 1 1') + CALL lmp%fix_external_set_vector_length('ext2', vec_length) +END SUBROUTINE f_lammps_setup_fix_external_callback + +SUBROUTINE f_lammps_setup_fix_external_array() BIND(C) + USE LIBLAMMPS + USE keepstuff, ONLY : lmp, demo_input, cont_input, pair_input + USE ext_stuff, ONLY : f3, f4 + IMPLICIT NONE + + CALL lmp%commands_list(demo_input) + CALL lmp%commands_list(cont_input) + CALL lmp%commands_list(pair_input) + CALL lmp%command('neigh_modify exclude group all all') + CALL lmp%command('fix ext3 all external pf/array 1') + CALL lmp%command('fix ext4 all external pf/array 1') + CALL lmp%command('thermo_style custom step pxx pe etotal') + CALL lmp%command('thermo_modify norm no') + CALL lmp%command('thermo 100') + f3 = lmp%fix_external_get_force('ext3') + f4 = lmp%fix_external_get_force('ext4') +END SUBROUTINE f_lammps_setup_fix_external_array + +SUBROUTINE f_lammps_set_fix_external_callbacks() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + USE ext_stuff + IMPLICIT NONE + INTEGER :: size_bigint, size_tagint, nlocal + + nlocal = lmp%extract_setting('nlocal') + + size_bigint = lmp%extract_setting('bigint') + size_tagint = lmp%extract_setting('tagint') + IF (size_bigint == 4_c_int .AND. size_tagint == 4_c_int) THEN + CALL lmp%set_fix_external_callback('ext1', f_callback_ss, lmp) + CALL lmp%set_fix_external_callback('ext2', f_callback2_ss, direction) + ELSE IF (size_bigint == 8_c_int .AND. size_tagint == 8_c_int) THEN + CALL lmp%set_fix_external_callback('ext1', f_callback_bb, lmp) + CALL lmp%set_fix_external_callback('ext2', f_callback2_bb, direction) + ELSE + CALL lmp%set_fix_external_callback('ext1', f_callback_sb, lmp) + CALL lmp%set_fix_external_callback('ext2', f_callback2_sb, direction) + END IF +END SUBROUTINE f_lammps_set_fix_external_callbacks + +SUBROUTINE f_lammps_get_force (i, ptr) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double, c_ptr, C_F_POINTER + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + INTEGER(c_int), INTENT(IN), VALUE :: i + TYPE(c_ptr), INTENT(IN), VALUE :: ptr + REAL(c_double), DIMENSION(:,:), POINTER :: force => NULL() + REAL(c_double), DIMENSION(:), POINTER :: f => NULL() + + CALL C_F_POINTER(ptr, f, [3]) + force = lmp%extract_atom('f') + f = force(:,i) +END SUBROUTINE f_lammps_get_force + +SUBROUTINE f_lammps_find_forces() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double, c_int, c_int64_t + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + USE ext_stuff, ONLY : f3, f4 + IMPLICIT NONE + INTEGER(c_int) :: size_tagint + INTEGER(c_int), DIMENSION(:), POINTER :: id + INTEGER(c_int64_t), DIMENSION(:), POINTER :: tag + + f3(:,:) = 0.0_c_double + f4(:,:) = 0.0_c_double + size_tagint = lmp%extract_setting('tagint') + IF (size_tagint == 4_c_int) THEN + id = lmp%extract_atom('id') + WHERE (id == 1_c_int) + f3(1,:) = 4.0_c_double + f3(2,:) = -4.0_c_double + f3(3,:) = 6.0_c_double + f4(1,:) = 10.0_c_double + f4(2,:) = -10.0_c_double + f4(3,:) = 12.0_c_double + ELSEWHERE + f3(1,:) = 5.0_c_double + f3(2,:) = -5.0_c_double + f3(3,:) = 7.0_c_double + f4(1,:) = 11.0_c_double + f4(2,:) = -11.0_c_double + f4(3,:) = 13.0_c_double + END WHERE + ELSE + tag = lmp%extract_atom('id') + WHERE (tag == 1_c_int64_t) + f3(1,:) = 4.0_c_double + f3(2,:) = -4.0_c_double + f3(3,:) = 6.0_c_double + f4(1,:) = 10.0_c_double + f4(2,:) = -10.0_c_double + f4(3,:) = 12.0_c_double + ELSEWHERE + f3(1,:) = 5.0_c_double + f3(2,:) = -5.0_c_double + f3(3,:) = 7.0_c_double + f4(1,:) = 11.0_c_double + f4(2,:) = -11.0_c_double + f4(3,:) = 13.0_c_double + END WHERE + END IF +END SUBROUTINE f_lammps_find_forces + +SUBROUTINE f_lammps_add_energy() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + + CALL lmp%fix_external_set_energy_global('ext3', -20.2_c_double); +END SUBROUTINE f_lammps_add_energy + +SUBROUTINE f_lammps_set_virial() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + + CALL lmp%fix_external_set_virial_global('ext4', [1.0_c_double, & + 2.0_c_double, 2.5_c_double, -1.0_c_double, -2.25_c_double, -3.02_c_double]) +END SUBROUTINE f_lammps_set_virial + +FUNCTION f_lammps_find_peratom_energy(i) RESULT(energy) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + INTEGER(c_int), INTENT(IN), VALUE :: i + REAL(c_double) :: energy + REAL(c_double), DIMENSION(:), POINTER :: e + + e = lmp%extract_compute('peratom', lmp%style%atom, lmp%type%vector) + energy = e(i) +END FUNCTION f_lammps_Find_peratom_energy + +SUBROUTINE f_lammps_find_peratom_virial(v, i) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + REAL(c_double), DIMENSION(6) :: v + INTEGER(c_int), INTENT(IN), VALUE :: i + REAL(c_double), DIMENSION(:,:), POINTER :: virial + + virial = lmp%extract_compute('vperatom', lmp%style%atom, lmp%type%array) + v = virial(:,i) +END SUBROUTINE f_lammps_find_peratom_virial + +SUBROUTINE f_lammps_fixexternal_set_vector() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + USE ext_stuff, ONLY : vec_length + IMPLICIT NONE + REAL(c_double), DIMENSION(vec_length) :: v + INTEGER :: i + DO i = 1, vec_length + v(i) = REAL(i, c_double) + CALL lmp%fix_external_set_vector('ext2', i, v(i)) + END DO +END SUBROUTINE f_lammps_fixexternal_set_vector diff --git a/unittest/fortran/test_fortran_gather_scatter.f90 b/unittest/fortran/test_fortran_gather_scatter.f90 index 2115a10394..6df1421bac 100644 --- a/unittest/fortran/test_fortran_gather_scatter.f90 +++ b/unittest/fortran/test_fortran_gather_scatter.f90 @@ -24,13 +24,16 @@ END SUBROUTINE f_lammps_close SUBROUTINE f_lammps_setup_gather_scatter() BIND(C) USE LIBLAMMPS - USE keepstuff, ONLY : lmp, big_input, cont_input, more_input + USE keepstuff, ONLY : lmp, big_input, cont_input, more_input, pair_input IMPLICIT NONE CALL lmp%command('atom_modify map array') CALL lmp%commands_list(big_input) CALL lmp%commands_list(cont_input) CALL lmp%commands_list(more_input) + CALL lmp%commands_list(pair_input) + CALL lmp%command('mass 1 1.0') + CALL lmp%command("compute pe all pe/atom") END SUBROUTINE f_lammps_setup_gather_scatter FUNCTION f_lammps_gather_atoms_mask(i) BIND(C) @@ -262,3 +265,90 @@ FUNCTION f_lammps_test_gather_bonds_big() BIND(C) RESULT(success) success = 0_c_int END IF END FUNCTION f_lammps_test_gather_bonds_big + +FUNCTION f_lammps_gather_pe_atom(i) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + INTEGER(c_int), INTENT(IN), VALUE :: i + REAL(c_double) :: f_lammps_gather_pe_atom + REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom + + CALL lmp%gather('c_pe', 1_c_int, pe_atom) + f_lammps_gather_pe_atom = pe_atom(i) +END FUNCTION f_lammps_gather_pe_atom + +FUNCTION f_lammps_gather_pe_atom_concat(i) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + INTEGER(c_int), INTENT(IN), VALUE :: i + REAL(c_double) :: f_lammps_gather_pe_atom_concat + REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom + INTEGER(c_int), DIMENSION(:), ALLOCATABLE :: tag + INTEGER :: j + + CALL lmp%gather_concat('id', 1_c_int, tag) + CALL lmp%gather_concat('c_pe', 1_c_int, pe_atom) + DO j = 1, SIZE(tag) + IF (tag(j) == i) THEN + f_lammps_gather_pe_atom_concat = pe_atom(j) + EXIT + END IF + END DO + f_lammps_gather_pe_atom_concat = pe_atom(i) +END FUNCTION f_lammps_gather_pe_atom_concat + +SUBROUTINE f_lammps_gather_pe_atom_subset(ids, pe) BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + INTEGER(c_int), INTENT(IN) :: ids(2) + REAL(c_double), INTENT(OUT) :: pe(2) + REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom + INTEGER(c_int) :: natoms + + natoms = NINT(lmp%get_natoms(), c_int) + CALL lmp%gather_subset('c_pe', 1, ids, pe_atom) + pe(1:natoms) = pe_atom +END SUBROUTINE f_lammps_gather_pe_atom_subset + +SUBROUTINE f_lammps_scatter_compute() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom + REAL(c_double) :: swap + + CALL lmp%gather('c_pe', 1_c_int, pe_atom) + + ! swap the computed energy of atoms 1 and 3 + swap = pe_atom(1) + pe_atom(1) = pe_atom(3) + pe_atom(3) = swap + + CALL lmp%scatter('c_pe', pe_atom) ! push the swap back to LAMMPS +END SUBROUTINE f_lammps_scatter_compute + +SUBROUTINE f_lammps_scatter_subset_compute() BIND(C) + USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double + USE LIBLAMMPS + USE keepstuff, ONLY : lmp + IMPLICIT NONE + INTEGER(c_int), PARAMETER :: ids(2) = [3,1] + REAL(c_double), DIMENSION(:), ALLOCATABLE :: pe_atom + REAL(c_double) :: swap + + CALL lmp%gather_subset('c_pe', 1_c_int, ids, pe_atom) + + ! swap the computed energy of atoms 1 and 3 + swap = pe_atom(1) + pe_atom(1) = pe_atom(2) + pe_atom(2) = swap + + CALL lmp%scatter_subset('c_pe', ids, pe_atom) ! push the swap back to LAMMPS +END SUBROUTINE f_lammps_scatter_subset_compute diff --git a/unittest/fortran/wrap_fixexternal.cpp b/unittest/fortran/wrap_fixexternal.cpp new file mode 100644 index 0000000000..694aceb009 --- /dev/null +++ b/unittest/fortran/wrap_fixexternal.cpp @@ -0,0 +1,194 @@ + +// unit tests for gathering and scattering 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_fix_external_callback(); +void f_lammps_setup_fix_external_array(); +void f_lammps_set_fix_external_callbacks(); +void f_lammps_get_force(int, double*); +void f_lammps_reverse_direction(); +void f_lammps_find_forces(); +void f_lammps_add_energy(); +void f_lammps_set_virial(); +double f_lammps_find_peratom_energy(int); +void f_lammps_find_peratom_virial(double[6], int); +void f_lammps_fixexternal_set_vector(); +} + +using namespace LAMMPS_NS; + +class LAMMPS_fixexternal : public ::testing::Test { +protected: + LAMMPS_NS::LAMMPS *lmp; + LAMMPS_fixexternal() = default; + ~LAMMPS_fixexternal() 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_fixexternal, callback) +{ + f_lammps_setup_fix_external_callback(); + f_lammps_set_fix_external_callbacks(); + lammps_command(lmp, "run 0"); + double f[3]; + f_lammps_get_force(1,f); + EXPECT_DOUBLE_EQ(f[0], 3.0); + EXPECT_DOUBLE_EQ(f[1], -3.0); + EXPECT_DOUBLE_EQ(f[2], 3.75); + f_lammps_get_force(2,f); + EXPECT_DOUBLE_EQ(f[0], -3.0); + EXPECT_DOUBLE_EQ(f[1], 3.0); + EXPECT_DOUBLE_EQ(f[2], -3.75); + + f_lammps_reverse_direction(); + f_lammps_set_fix_external_callbacks(); + lammps_command(lmp, "run 0"); + f_lammps_get_force(1,f); + EXPECT_DOUBLE_EQ(f[0], -1.0); + EXPECT_DOUBLE_EQ(f[1], 1.0); + EXPECT_DOUBLE_EQ(f[2], -1.25); + f_lammps_get_force(2,f); + EXPECT_DOUBLE_EQ(f[0], 1.0); + EXPECT_DOUBLE_EQ(f[1], -1.0); + EXPECT_DOUBLE_EQ(f[2], 1.25); +}; + +TEST_F(LAMMPS_fixexternal, array) +{ + f_lammps_setup_fix_external_array(); + double **f; + f = (double**) lammps_extract_atom(lmp, "f"); + f_lammps_find_forces(); + lammps_command(lmp, "run 0"); + EXPECT_DOUBLE_EQ(f[0][0], 14.0); + EXPECT_DOUBLE_EQ(f[0][1], -14.0); + EXPECT_DOUBLE_EQ(f[0][2], 18.0); + EXPECT_DOUBLE_EQ(f[1][0], 16.0); + EXPECT_DOUBLE_EQ(f[1][1], -16.0); + EXPECT_DOUBLE_EQ(f[1][2], 20.0); +}; + +TEST_F(LAMMPS_fixexternal, energy_global) +{ + f_lammps_setup_fix_external_array(); + double energy; + f_lammps_add_energy(); + lammps_command(lmp, "run 0"); + energy = lammps_get_thermo(lmp, "etotal"); + EXPECT_DOUBLE_EQ(energy, -20.2); +}; + +TEST_F(LAMMPS_fixexternal, virial_global) +{ + f_lammps_setup_fix_external_array(); + double virial[6], volume; + f_lammps_set_virial(); + lammps_command(lmp, "run 0"); + volume = lammps_get_thermo(lmp, "vol"); + virial[0] = lammps_get_thermo(lmp, "pxx"); + virial[1] = lammps_get_thermo(lmp, "pyy"); + virial[2] = lammps_get_thermo(lmp, "pzz"); + virial[3] = lammps_get_thermo(lmp, "pxy"); + virial[4] = lammps_get_thermo(lmp, "pxz"); + virial[5] = lammps_get_thermo(lmp, "pyz"); + EXPECT_DOUBLE_EQ(virial[0], 1.0/volume); + EXPECT_DOUBLE_EQ(virial[1], 2.0/volume); + EXPECT_DOUBLE_EQ(virial[2], 2.5/volume); + EXPECT_DOUBLE_EQ(virial[3], -1.0/volume); + EXPECT_DOUBLE_EQ(virial[4], -2.25/volume); + EXPECT_DOUBLE_EQ(virial[5], -3.02/volume); +}; + +TEST_F(LAMMPS_fixexternal, energy_peratom) +{ + f_lammps_setup_fix_external_callback(); + f_lammps_set_fix_external_callbacks(); + lammps_command(lmp, "compute peratom all pe/atom"); + double energy; + lammps_command(lmp, "run 0"); + int nlocal = lammps_extract_setting(lmp, "nlocal"); + for (int i = 1; i <= nlocal; i++) + { + energy = f_lammps_find_peratom_energy(i); + if (i == 1) + EXPECT_DOUBLE_EQ(energy, 1.0); + else + EXPECT_DOUBLE_EQ(energy, 10.0); + } +}; + +TEST_F(LAMMPS_fixexternal, virial_peratom) +{ + f_lammps_setup_fix_external_callback(); + f_lammps_set_fix_external_callbacks(); + lammps_command(lmp, "compute vperatom all stress/atom NULL"); + double virial[6]; + lammps_command(lmp, "run 0"); + int nlocal = lammps_extract_setting(lmp, "nlocal"); + for (int i = 1; i <= nlocal; i++) + { + f_lammps_find_peratom_virial(virial, i); + if (i == 1) + { + EXPECT_DOUBLE_EQ(virial[0], -1.0); + EXPECT_DOUBLE_EQ(virial[1], -2.0); + EXPECT_DOUBLE_EQ(virial[2], 1.0); + EXPECT_DOUBLE_EQ(virial[3], 2.0); + EXPECT_DOUBLE_EQ(virial[4], -3.0); + EXPECT_DOUBLE_EQ(virial[5], 3.0); + } + else + { + EXPECT_DOUBLE_EQ(virial[0], -10.0); + EXPECT_DOUBLE_EQ(virial[1], -20.0); + EXPECT_DOUBLE_EQ(virial[2], 10.0); + EXPECT_DOUBLE_EQ(virial[3], 20.0); + EXPECT_DOUBLE_EQ(virial[4], -30.0); + EXPECT_DOUBLE_EQ(virial[5], 30.0); + } + } +}; + +TEST_F(LAMMPS_fixexternal, vector) +{ + f_lammps_setup_fix_external_callback(); + f_lammps_set_fix_external_callbacks(); + f_lammps_fixexternal_set_vector(); + lammps_command(lmp, "run 0"); + double *v; + for (int i = 0; i < 8; i++) + { + v = (double*) lammps_extract_fix(lmp, "ext2", LMP_STYLE_GLOBAL, + LMP_TYPE_VECTOR, i, 1); + EXPECT_DOUBLE_EQ(i+1, *v); + std::free(v); + } +}; diff --git a/unittest/fortran/wrap_gather_scatter.cpp b/unittest/fortran/wrap_gather_scatter.cpp index f200e91b48..5eecf763a6 100644 --- a/unittest/fortran/wrap_gather_scatter.cpp +++ b/unittest/fortran/wrap_gather_scatter.cpp @@ -3,6 +3,7 @@ #include "lammps.h" #include "library.h" +#include "atom.h" #include #include #include @@ -26,6 +27,11 @@ void f_lammps_scatter_atoms_positions(); void f_lammps_setup_gather_bonds(); int f_lammps_test_gather_bonds_small(); int f_lammps_test_gather_bonds_big(); +double f_lammps_gather_pe_atom(int); +double f_lammps_gather_pe_atom_concat(int); +void f_lammps_gather_pe_atom_subset(int*, double*); +void f_lammps_scatter_compute(); +void f_lammps_scatter_subset_compute(); } using namespace LAMMPS_NS; @@ -216,3 +222,113 @@ TEST_F(LAMMPS_gather_scatter, gather_bonds) EXPECT_EQ(f_lammps_test_gather_bonds_small(), 1); #endif }; + +TEST_F(LAMMPS_gather_scatter, gather_compute) +{ +#ifdef LAMMPS_BIGBIG + GTEST_SKIP(); +#else + f_lammps_setup_gather_scatter(); + lammps_command(lmp, "run 0"); + int natoms = lmp->atom->natoms; + int *tag = lmp->atom->tag; + double *pe = (double*) lammps_extract_compute(lmp, "pe", LMP_STYLE_ATOM, + LMP_TYPE_VECTOR); + for (int i = 0; i < natoms; i++) + EXPECT_DOUBLE_EQ(f_lammps_gather_pe_atom(tag[i]), pe[i]); +#endif +}; + +TEST_F(LAMMPS_gather_scatter, gather_compute_concat) +{ +#ifdef LAMMPS_BIGBIG + GTEST_SKIP(); +#else + f_lammps_setup_gather_scatter(); + lammps_command(lmp, "run 0"); + int natoms = lmp->atom->natoms; + int *tag = lmp->atom->tag; + double *pe = (double*) lammps_extract_compute(lmp, "pe", LMP_STYLE_ATOM, + LMP_TYPE_VECTOR); + for (int i = 0; i < natoms; i++) + EXPECT_DOUBLE_EQ(f_lammps_gather_pe_atom(tag[i]), pe[i]); +#endif +}; + +TEST_F(LAMMPS_gather_scatter, gather_compute_subset) +{ +#ifdef LAMMPS_BIGBIG + GTEST_SKIP(); +#else + f_lammps_setup_gather_scatter(); + lammps_command(lmp, "run 0"); + int ids[2] = {3, 1}; + int *tag = lmp->atom->tag; + double pe[2] = {0.0, 0.0}; + int nlocal = lammps_extract_setting(lmp, "nlocal"); + double *pa_pe = (double*) lammps_extract_compute(lmp, "pe", LMP_STYLE_ATOM, + LMP_TYPE_VECTOR); + + for (int i = 0; i < nlocal; i++) { + if(tag[i] == ids[0]) pe[0] = pa_pe[i]; + if(tag[i] == ids[1]) pe[1] = pa_pe[i]; + } + + double ftn_pe[2]; + f_lammps_gather_pe_atom_subset(ids, ftn_pe); + EXPECT_DOUBLE_EQ(ftn_pe[0], pe[0]); + EXPECT_DOUBLE_EQ(ftn_pe[1], pe[1]); +#endif +}; + +TEST_F(LAMMPS_gather_scatter, scatter_compute) +{ +#ifdef LAMMPS_BIGBIG + GTEST_SKIP(); +#else + f_lammps_setup_gather_scatter(); + int natoms = lmp->atom->natoms; + double *pe = new double[natoms]; + lammps_command(lmp, "run 0"); + lammps_gather(lmp, "c_pe", 1, 1, pe); + double *old_pe = new double[natoms]; + for (int i = 0; i < natoms; i++) + old_pe[i] = pe[i]; + EXPECT_DOUBLE_EQ(pe[0], old_pe[0]); + EXPECT_DOUBLE_EQ(pe[1], old_pe[1]); + EXPECT_DOUBLE_EQ(pe[2], old_pe[2]); + f_lammps_scatter_compute(); + lammps_gather(lmp, "c_pe", 1, 1, pe); + EXPECT_DOUBLE_EQ(pe[0], old_pe[2]); + EXPECT_DOUBLE_EQ(pe[1], old_pe[1]); + EXPECT_DOUBLE_EQ(pe[2], old_pe[0]); + delete[] old_pe; + delete[] pe; +#endif +}; + +TEST_F(LAMMPS_gather_scatter, scatter_subset_compute) +{ +#ifdef LAMMPS_BIGBIG + GTEST_SKIP(); +#else + f_lammps_setup_gather_scatter(); + int natoms = lmp->atom->natoms; + double *pe = new double[natoms]; + lammps_command(lmp, "run 0"); + lammps_gather(lmp, "c_pe", 1, 1, pe); + double *old_pe = new double[natoms]; + for (int i = 0; i < natoms; i++) + old_pe[i] = pe[i]; + EXPECT_DOUBLE_EQ(pe[0], old_pe[0]); + EXPECT_DOUBLE_EQ(pe[1], old_pe[1]); + EXPECT_DOUBLE_EQ(pe[2], old_pe[2]); + f_lammps_scatter_subset_compute(); + lammps_gather(lmp, "c_pe", 1, 1, pe); + EXPECT_DOUBLE_EQ(pe[0], old_pe[2]); + EXPECT_DOUBLE_EQ(pe[1], old_pe[1]); + EXPECT_DOUBLE_EQ(pe[2], old_pe[0]); + delete[] old_pe; + delete[] pe; +#endif +};