diff --git a/doc/src/Fortran.rst b/doc/src/Fortran.rst index 874b57cbd7..dc5001f4c9 100644 --- a/doc/src/Fortran.rst +++ b/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 ` 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 ` 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 - part of the callback function. For this to work, the handle to the LAMMPS - object must be passed as the *ptr* argument when registering the callback - function. + 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:: @@ -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` diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index c5ccc17b4e..fc318091bf 100644 --- a/fortran/lammps.f90 +++ b/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) diff --git a/unittest/fortran/test_fortran_fixexternal.f90 b/unittest/fortran/test_fortran_fixexternal.f90 index 1da8909627..ca08e3c4d9 100644 --- a/unittest/fortran/test_fortran_fixexternal.f90 +++ b/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 diff --git a/unittest/fortran/wrap_fixexternal.cpp b/unittest/fortran/wrap_fixexternal.cpp index 282b95a299..aa6cb2b921 100644 --- a/unittest/fortran/wrap_fixexternal.cpp +++ b/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); +};