Fixed bug and wrote unit tests for fix_external_array functions

doc/src/Fortran.rst

@@ -2254,7 +2254,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type
    energy units as determined by the current :doc:`units <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 scalar property of the fix external
+   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 <fix_external>` for more
@@ -2307,10 +2307,11 @@ Procedures Bound to the :f:type:`lammps` Derived Type
    .. 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 code as
+   the per-atom energy contribution due to the fix from the external program as
-   part of the callback function.  For this to work, the handle to the LAMMPS
+   part of the callback function.  For this to work, the LAMMPS object must be
-   object must be passed as the *ptr* argument when registering the callback
+   passed as part of the *caller* argument when registering the callback
-   function.
+   function, or the callback function must otherwise have access to the
+   LAMMPS object, such as through a module-based pointer.
 
    .. note::
 
@@ -2325,7 +2326,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type
    external code.
 
    :p character(len=*) id: fix ID of the fix external instance
-   :p real(c_double) eng [dimension(:)]: array of length nlocal containing
+   :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`
 

fortran/lammps.f90

@@ -2618,9 +2618,9 @@ CONTAINS
 
   ! equivalent function to lammps_fix_external_get_force
   FUNCTION lmp_fix_external_get_force(self, id) RESULT(fexternal)
-    CLASS(lammps), INTENT(IN) :: self
+    CLASS(lammps), TARGET, INTENT(IN) :: self
     CHARACTER(LEN=*), INTENT(IN) :: id
-    REAL(c_double), DIMENSION(:,:), POINTER :: fexternal
+    TYPE(lammps_fix_data) :: fexternal
     TYPE(c_ptr) :: ptr, Cid
     TYPE(c_ptr), DIMENSION(:), POINTER :: f
     INTEGER(c_int) :: nmax
@@ -2629,14 +2629,16 @@ CONTAINS
     ptr = lammps_fix_external_get_force(self%handle, Cid)
     nmax = lmp_extract_setting(self, 'nmax')
     CALL C_F_POINTER(ptr, f, [nmax])
-    CALL C_F_POINTER(f(1), fexternal, [3, 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(OUT) :: eng
+    REAL(c_double), INTENT(IN) :: eng
     TYPE(c_ptr) :: Cid
 
     Cid = f2c_string(id)

unittest/fortran/test_fortran_fixexternal.f90

@@ -5,6 +5,7 @@ MODULE ext_stuff
   IMPLICIT NONE
 
   REAL(c_double), SAVE :: direction = 1.0_c_double
+  REAL(c_double), DIMENSION(:,:), POINTER, SAVE :: f3 => NULL(), f4 => NULL()
 
 CONTAINS
 
@@ -18,16 +19,26 @@ CONTAINS
     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
 
     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
     ELSEWHERE
       f(1,:) = -1.0_c_double
       f(2,:) = +1.0_c_double
       f(3,:) = -1.25_c_double
+      e = 10.0_c_double
     END WHERE
+    SELECT TYPE (instance)
+      CLASS IS (lammps)
+        CALL instance%fix_external_set_energy_peratom('ext1', e)
+      CLASS DEFAULT
+        WRITE(0,*) 'UMM...this should never happen.'
+        STOP 1
+    END SELECT
   END SUBROUTINE f_callback_ss
 
   SUBROUTINE f_callback_sb(instance, timestep, id, x, f)
@@ -36,16 +47,26 @@ CONTAINS
     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
 
     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
     ELSEWHERE
       f(1,:) = -1.0_c_double
       f(2,:) = +1.0_c_double
       f(3,:) = -1.25_c_double
+      e = 10.0_c_double
     END WHERE
+    SELECT TYPE (instance)
+      CLASS IS (lammps)
+        CALL instance%fix_external_set_energy_peratom('ext1', e)
+      CLASS DEFAULT
+        WRITE(0,*) 'UMM...this should never happen.'
+        STOP 1
+    END SELECT
   END SUBROUTINE f_callback_sb
 
   SUBROUTINE f_callback_bb(instance, timestep, id, x, f)
@@ -54,16 +75,26 @@ CONTAINS
     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
 
     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
     ELSEWHERE
       f(1,:) = -1.0_c_double
       f(2,:) = +1.0_c_double
       f(3,:) = -1.25_c_double
+      e = 10.0_c_double
     END WHERE
+    SELECT TYPE (instance)
+      CLASS IS (lammps)
+        CALL instance%fix_external_set_energy_peratom('ext1', e)
+      CLASS DEFAULT
+        WRITE(0,*) 'UMM...this should never happen.'
+        STOP 1
+    END SELECT
   END SUBROUTINE f_callback_bb
 
   SUBROUTINE f_callback2_ss(entity, timestep, id, x, f)
@@ -158,7 +189,7 @@ FUNCTION f_lammps_with_args() BIND(C)
 END FUNCTION f_lammps_with_args
 
 SUBROUTINE f_lammps_close() BIND(C)
-  USE ISO_C_BINDING, ONLY: c_null_ptr
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY: c_null_ptr
   USE liblammps
   USE keepstuff, ONLY: lmp
   IMPLICIT NONE
@@ -167,7 +198,7 @@ SUBROUTINE f_lammps_close() BIND(C)
   lmp%handle = c_null_ptr
 END SUBROUTINE f_lammps_close
 
-SUBROUTINE f_lammps_setup_fix_external() BIND(C)
+SUBROUTINE f_lammps_setup_fix_external_callback() BIND(C)
   USE LIBLAMMPS
   USE keepstuff, ONLY : lmp, demo_input, cont_input, pair_input
   IMPLICIT NONE
@@ -178,10 +209,29 @@ SUBROUTINE f_lammps_setup_fix_external() BIND(C)
   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')
-END SUBROUTINE f_lammps_setup_fix_external
+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 ISO_C_BINDING, ONLY : c_int
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int
   USE LIBLAMMPS
   USE keepstuff, ONLY : lmp
   USE ext_stuff
@@ -193,19 +243,19 @@ SUBROUTINE f_lammps_set_fix_external_callbacks() BIND(C)
   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)
+    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)
+    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)
+    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 ISO_C_BINDING, ONLY : c_int, c_double, c_ptr, C_F_POINTER
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : c_int, c_double, c_ptr, C_F_POINTER
   USE LIBLAMMPS
   USE keepstuff, ONLY : lmp
   IMPLICIT NONE
@@ -218,3 +268,69 @@ SUBROUTINE f_lammps_get_force (i, ptr) BIND(C)
   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
+    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

unittest/fortran/wrap_fixexternal.cpp

@@ -15,10 +15,14 @@
 extern "C" {
 void *f_lammps_with_args();
 void f_lammps_close();
-void f_lammps_setup_fix_external();
+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();
 }
 
 using namespace LAMMPS_NS;
@@ -48,7 +52,7 @@ protected:
 
 TEST_F(LAMMPS_fixexternal, callback)
 {
-    f_lammps_setup_fix_external();
+    f_lammps_setup_fix_external_callback();
     f_lammps_set_fix_external_callbacks();
     lammps_command(lmp, "run 0");
     double f[3];
@@ -73,3 +77,61 @@ TEST_F(LAMMPS_fixexternal, callback)
     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();
+    double energy;
+    lammps_command(lmp, "run 0");
+/* FIXME: the per-atom energy is NOT summed up by this function! We need
+   another test. */
+    energy = lammps_get_thermo(lmp, "pe");
+    EXPECT_DOUBLE_EQ(energy, 11.0);
+};