diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index fc318091bf..18b8fa89a3 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -2682,8 +2682,8 @@ CONTAINS CHARACTER(LEN=*), INTENT(IN) :: id REAL(c_double), DIMENSION(:,:), TARGET, INTENT(IN) :: virial TYPE(c_ptr) :: Cid, Cvirial - TYPE(c_ptr), TARGET :: Cptr - INTEGER(c_int) :: nlocal + 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 @@ -2692,10 +2692,14 @@ CONTAINS &[Fortran/fix_external_set_energy_peratom]') END IF Cid = f2c_string(id) - Cptr = C_LOC(virial(1,1)) - Cvirial = C_LOC(Cptr) + 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) diff --git a/unittest/fortran/test_fortran_fixexternal.f90 b/unittest/fortran/test_fortran_fixexternal.f90 index ca08e3c4d9..61fea1be9d 100644 --- a/unittest/fortran/test_fortran_fixexternal.f90 +++ b/unittest/fortran/test_fortran_fixexternal.f90 @@ -4,6 +4,7 @@ MODULE ext_stuff 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() @@ -20,21 +21,35 @@ CONTAINS 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(0,*) 'UMM...this should never happen.' STOP 1 @@ -48,21 +63,35 @@ CONTAINS 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(0,*) 'UMM...this should never happen.' STOP 1 @@ -76,21 +105,35 @@ CONTAINS 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(0,*) 'UMM...this should never happen.' STOP 1 @@ -201,6 +244,7 @@ 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) @@ -209,6 +253,7 @@ SUBROUTINE f_lammps_setup_fix_external_callback() 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') + 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) @@ -305,6 +350,9 @@ SUBROUTINE f_lammps_find_forces() BIND(C) 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 @@ -334,3 +382,43 @@ SUBROUTINE f_lammps_set_virial() BIND(C) 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/wrap_fixexternal.cpp b/unittest/fortran/wrap_fixexternal.cpp index aa6cb2b921..694aceb009 100644 --- a/unittest/fortran/wrap_fixexternal.cpp +++ b/unittest/fortran/wrap_fixexternal.cpp @@ -23,6 +23,9 @@ 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; @@ -128,10 +131,64 @@ 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"); -/* 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); + 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); + } };