git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9003 f3b2605a-c512-4ea7-a41b-209d697bcdaa

examples/COUPLE/fortran2/LAMMPS.F90

@@ -57,7 +57,7 @@ module LAMMPS
       lammps_extract_atom, lammps_extract_compute, lammps_extract_fix, &
       lammps_extract_variable, lammps_get_natoms, lammps_gather_atoms, &
       lammps_scatter_atoms
-   public :: lammps_instance
+   public :: lammps_instance, C_ptr, C_double, C_int
 
    !! Functions supplemental to the prototypes in library.h. {{{1
    !! The function definitions (in C++) are contained in LAMMPS-wrapper.cpp.
@@ -224,52 +224,40 @@ module LAMMPS
 
    ! Generic functions for the wrappers below {{{1
 
-   ! Check the dimensions of the arrays these return; they are not always
-   ! easy to find.  Note that I consider returning pointers to arbitrary
-   ! memory locations with no information as to array size/shape to be
-   ! extremely sloppy and error-prone.  It would appear the Fortran standards
-   ! committee would agree, as they chose not to allow that sort of nonsense.
-
    interface lammps_extract_global
-      module procedure lammps_extract_global_i, lammps_extract_global_r, &
+      module procedure lammps_extract_global_i, &
          lammps_extract_global_dp
    end interface lammps_extract_global
 
    interface lammps_extract_atom
-      module procedure lammps_extract_atom_ia, lammps_extract_atom_ra, &
+      module procedure lammps_extract_atom_ia, &
-         lammps_extract_atom_dpa, lammps_extract_atom_dp2a, &
+         lammps_extract_atom_dpa, &
-         lammps_extract_atom_r2a
+         lammps_extract_atom_dp2a
    end interface lammps_extract_atom
 
    interface lammps_extract_compute
-      module procedure lammps_extract_compute_r, lammps_extract_compute_dp, &
+      module procedure lammps_extract_compute_dp, &
-         lammps_extract_compute_ra, lammps_extract_compute_dpa, &
+         lammps_extract_compute_dpa, &
-         lammps_extract_compute_r2a, lammps_extract_compute_dp2a
+         lammps_extract_compute_dp2a
    end interface lammps_extract_compute
 
    interface lammps_extract_fix
-      module procedure lammps_extract_fix_r, lammps_extract_fix_dp, &
+      module procedure lammps_extract_fix_dp, &
-         lammps_extract_fix_ra, lammps_extract_fix_dpa, &
+         lammps_extract_fix_dpa, &
-         lammps_extract_fix_r2a, lammps_extract_fix_dp2a
+         lammps_extract_fix_dp2a
    end interface lammps_extract_fix
 
    interface lammps_extract_variable
-      module procedure lammps_extract_variable_i, &
+      module procedure lammps_extract_variable_dp, &
-         lammps_extract_variable_dp, &
-         lammps_extract_variable_r, &
-         lammps_extract_variable_ra, &
-         lammps_extract_variable_ia, &
          lammps_extract_variable_dpa
    end interface lammps_extract_variable
 
    interface lammps_gather_atoms
-      module procedure lammps_gather_atoms_ia, lammps_gather_atoms_dpa, &
+      module procedure lammps_gather_atoms_ia, lammps_gather_atoms_dpa
-         lammps_gather_atoms_ra
    end interface lammps_gather_atoms
 
    interface lammps_scatter_atoms
-      module procedure lammps_scatter_atoms_ia, lammps_scatter_atoms_dpa, &
+      module procedure lammps_scatter_atoms_ia, lammps_scatter_atoms_dpa
-         lammps_scatter_atoms_ra
    end interface lammps_scatter_atoms
 
 contains !! Wrapper functions local to this module {{{1
@@ -336,38 +324,21 @@ contains !! Wrapper functions local to this module {{{1
       global = lammps_actual_extract_global (ptr, Cname)
    end function lammps_extract_global_Cptr
    subroutine lammps_extract_global_i (global, ptr, name)
-      integer, intent(out) :: global
+      integer (C_int), pointer, intent(out) :: global
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: name
       type (C_ptr) :: Cptr
-      integer (C_int), pointer :: Fptr
       Cptr = lammps_extract_global_Cptr (ptr, name)
-      call C_F_pointer (Cptr, Fptr)
+      call C_F_pointer (Cptr, global)
-      global = Fptr
-      nullify (Fptr)
    end subroutine lammps_extract_global_i
    subroutine lammps_extract_global_dp (global, ptr, name)
-      double precision, intent(out) :: global
+      real (C_double), pointer, intent(out) :: global
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: name
       type (C_ptr) :: Cptr
-      real (C_double), pointer :: Fptr
       Cptr = lammps_extract_global_Cptr (ptr, name)
-      call C_F_pointer (Cptr, Fptr)
+      call C_F_pointer (Cptr, global)
-      global = Fptr
-      nullify (Fptr)
    end subroutine lammps_extract_global_dp
-   subroutine lammps_extract_global_r (global, ptr, name)
-      real :: global
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      type (C_ptr) :: Cptr
-      real (C_double), pointer :: Fptr
-      Cptr = lammps_extract_global_Cptr (ptr, name)
-      call C_F_pointer (Cptr, Fptr)
-      global = real (Fptr)
-      nullify (Fptr)
-   end subroutine lammps_extract_global_r
 
 !-----------------------------------------------------------------------------
 
@@ -381,92 +352,69 @@ contains !! Wrapper functions local to this module {{{1
       atom = lammps_actual_extract_atom (ptr, Cname)
    end function lammps_extract_atom_Cptr
    subroutine lammps_extract_atom_ia (atom, ptr, name)
-      integer, dimension(:), allocatable, intent(out) :: atom
+      integer (C_int), dimension(:), pointer, intent(out) :: atom
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: name
       type (C_ptr) :: Cptr
-      integer (C_int), dimension(:), pointer :: Fptr
+      integer (C_int), pointer :: nelements
-      integer :: nelements
       call lammps_extract_global_i (nelements, ptr, 'nlocal')
       Cptr = lammps_extract_atom_Cptr (ptr, name)
-      call C_F_pointer (Cptr, Fptr, (/nelements/))
+      call C_F_pointer (Cptr, atom, (/nelements/))
-      if ( .not. associated (Fptr) ) return
-      allocate (atom(nelements))
-      atom = Fptr
-      nullify (Fptr)
    end subroutine lammps_extract_atom_ia
    subroutine lammps_extract_atom_dpa (atom, ptr, name)
-      double precision, dimension(:), allocatable, intent(out) :: atom
+      real (C_double), dimension(:), pointer, intent(out) :: atom
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: name
       type (C_ptr) :: Cptr
-      real (C_double), dimension(:), pointer :: Fptr
+      integer (C_int), pointer :: nlocal
       integer :: nelements
+      real (C_double), dimension(:), pointer :: Fptr
       if ( name == 'mass' ) then
-         nelements = lammps_get_ntypes (ptr)
+         nelements = lammps_get_ntypes (ptr) + 1
-      else if ( name == 'x' .or. name == 'v' .or. name == 'f' ) then
+      else if ( name == 'x' .or. name == 'v' .or. name == 'f' .or. &
-         ! We should not be getting 'x' or 'v' or 'f' here!
+                name == 'mu' .or. name == 'omega' .or. name == 'torque' .or. &
+                name == 'angmom' ) then
+         ! We should not be getting a rank-2 array here!
          call lammps_error_all (ptr, FLERR, 'You cannot extract those atom&
-            & data (x, v, or f) into a rank 1 array.')
+            & data (' // trim(name) // ') into a rank 1 array.')
          return
       else
          ! Everything else we can get is probably nlocal units long
-         call lammps_extract_global_i (nelements, ptr, 'nlocal')
+         call lammps_extract_global_i (nlocal, ptr, 'nlocal')
+         nelements = nlocal
       end if
       Cptr = lammps_extract_atom_Cptr (ptr, name)
-      if ( name == 'mass' ) then
-         call C_F_pointer (Cptr, Fptr, (/nelements + 1/))
-         if ( .not. associated (Fptr) ) return
-         allocate (atom(nelements))
-         atom = Fptr(2:) ! LAMMPS starts numbering at 1 (C does not)
-      else
       call C_F_pointer (Cptr, Fptr, (/nelements/))
-         if ( .not. associated (Fptr) ) return
+      if ( name == 'mass' ) then
-         allocate (atom(nelements))
+         atom(0:) => Fptr
-         atom = Fptr
+      else
+         atom => Fptr
       end if
-      nullify (Fptr)
    end subroutine lammps_extract_atom_dpa
-   subroutine lammps_extract_atom_ra (atom, ptr, name)
-      real, dimension(:), allocatable, intent(out) :: atom
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      double precision, dimension(:), allocatable :: d_atom
-      call lammps_extract_atom_dpa (d_atom, ptr, name)
-      allocate (atom(size(d_atom)))
-      atom = real(d_atom)
-      deallocate (d_atom)
-   end subroutine lammps_extract_atom_ra
    subroutine lammps_extract_atom_dp2a (atom, ptr, name)
-      double precision, dimension(:,:), allocatable, intent(out) :: atom
+      real (C_double), dimension(:,:), pointer, intent(out) :: atom
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: name
       type (C_ptr) :: Cptr
-      integer :: nelements
+      type (C_ptr), pointer, dimension(:) :: Catom
-      if ( name /= 'x' .and. name /= 'v' .and. name /= 'f' ) then
+      integer (C_int), pointer :: nelements
-         call lammps_error_all (ptr, FLERR, 'You cannot extract ' // name // &
+      if ( name /= 'x' .and. name /= 'v' .and. name /= 'f' .and. &
-            ' into a rank 2 array.')
+           name /= 'mu' .and. name /= 'omega' .and. name /= 'tandque' .and. &
+           name /= 'angmom' ) then
+         ! We should not be getting a rank-2 array here!
+         call lammps_error_all (ptr, FLERR, 'You cannot extract those atom&
+            & data (' // trim(name) // ') into a rank 2 array.')
          return
       end if
       Cptr = lammps_extract_atom_Cptr (ptr, name)
       call lammps_extract_global_i (nelements, ptr, 'nlocal')
-      allocate (atom(nelements,3))
+      ! Catom will now be the array of void* pointers that the void** pointer
-      atom = Cdoublestar_to_2darray (Cptr, nelements, 3)
+      ! pointed to.  Catom(1) is now the pointer to the first element.
+      call C_F_pointer (Cptr, Catom, (/nelements/))
+      ! Now get the actual array, which has its shape transposed from what we
+      ! might think of it in C
+      call C_F_pointer (Catom(1), atom, (/3, nelements/))
    end subroutine lammps_extract_atom_dp2a
-   subroutine lammps_extract_atom_r2a (atom, ptr, name)
-      real, dimension(:,:), allocatable, intent(out) :: atom
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      double precision, dimension(:,:), allocatable :: d_atom
-      call lammps_extract_atom_dp2a (d_atom, ptr, name)
-      if ( allocated (d_atom) ) then
-         allocate (atom(size(d_atom,1), size(d_atom,2)))
-      else
-         return
-      end if
-      atom = real(d_atom)
-      deallocate (d_atom)
-   end subroutine lammps_extract_atom_r2a
 
 !-----------------------------------------------------------------------------
 
@@ -484,12 +432,11 @@ contains !! Wrapper functions local to this module {{{1
       compute = lammps_actual_extract_compute (ptr, Cid, Cstyle, Ctype)
    end function lammps_extract_compute_Cptr
    subroutine lammps_extract_compute_dp (compute, ptr, id, style, type)
-      double precision, intent(out) :: compute
+      real (C_double), pointer, intent(out) :: compute
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: id
       integer, intent(in) :: style, type
       type (C_ptr) :: Cptr
-      real (C_double), pointer :: Fptr
       ! The only valid values of (style,type) are (0,0) for scalar 'compute'
       if ( style /= 0 ) then
          call lammps_error_all (ptr, FLERR, 'You cannot pack per-atom/local&
@@ -506,27 +453,14 @@ contains !! Wrapper functions local to this module {{{1
          return
       end if
       Cptr = lammps_extract_compute_Cptr (ptr, id, style, type)
-      call C_F_pointer (Cptr, Fptr)
+      call C_F_pointer (Cptr, compute)
-      compute = Fptr
-      nullify (Fptr)
-      ! C pointer should not be freed!
    end subroutine lammps_extract_compute_dp
-   subroutine lammps_extract_compute_r (compute, ptr, id, style, type)
-      real, intent(out) :: compute
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: id
-      integer, intent(in) :: style, type
-      double precision :: d_compute
-      call lammps_extract_compute_dp (d_compute, ptr, id, style, type)
-      compute = real(d_compute)
-   end subroutine lammps_extract_compute_r
    subroutine lammps_extract_compute_dpa (compute, ptr, id, style, type)
-      double precision, dimension(:), allocatable, intent(out) :: compute
+      real (C_double), dimension(:), pointer, intent(out) :: compute
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: id
       integer, intent(in) :: style, type
       type (C_ptr) :: Cptr
-      real (C_double), dimension(:), pointer :: Fptr
       integer :: nelements
       ! Check for the correct dimensionality
       if ( type == 0 ) then
@@ -539,30 +473,16 @@ contains !! Wrapper functions local to this module {{{1
          return
       end if
       nelements = lammps_extract_compute_vectorsize (ptr, id, style)
-      allocate (compute(nelements))
       Cptr = lammps_extract_compute_Cptr (ptr, id, style, type)
-      call C_F_pointer (Cptr, Fptr, (/nelements/))
+      call C_F_pointer (Cptr, compute, (/nelements/))
-      compute = Fptr
-      nullify (Fptr)
-      ! C pointer should not be freed
    end subroutine lammps_extract_compute_dpa
-   subroutine lammps_extract_compute_ra (compute, ptr, id, style, type)
-      real, dimension(:), allocatable, intent(out) :: compute
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: id
-      integer, intent(in) :: style, type
-      double precision, dimension(:), allocatable :: d_compute
-      call lammps_extract_compute_dpa (d_compute, ptr, id, style, type)
-      allocate (compute(size(d_compute)))
-      compute = real(d_compute)
-      deallocate (d_compute)
-   end subroutine lammps_extract_compute_ra
    subroutine lammps_extract_compute_dp2a (compute, ptr, id, style, type)
-      double precision, dimension(:,:), allocatable, intent(out) :: compute
+      real (C_double), dimension(:,:), pointer, intent(out) :: compute
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: id
       integer, intent(in) :: style, type
       type (C_ptr) :: Cptr
+      type (C_ptr), pointer, dimension(:) :: Ccompute
       integer :: nr, nc
       ! Check for the correct dimensionality
       if ( type == 0 ) then
@@ -575,22 +495,10 @@ contains !! Wrapper functions local to this module {{{1
          return
       end if
       call lammps_extract_compute_arraysize (ptr, id, style, nr, nc)
-      allocate (compute(nr, nc))
+      call C_F_pointer (Cptr, Ccompute, (/nr/))
-      Cptr = lammps_extract_compute_Cptr (ptr, id, style, type)
+      ! Note that the matrix is transposed, from Fortran's perspective
-      compute = Cdoublestar_to_2darray (Cptr, nr, nc)
+      call C_F_pointer (Ccompute(1), compute, (/nc, nr/))
-      ! C pointer should not be freed
    end subroutine lammps_extract_compute_dp2a
-   subroutine lammps_extract_compute_r2a (compute, ptr, id, style, type)
-      real, dimension(:,:), allocatable, intent(out) :: compute
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: id
-      integer, intent(in) :: style, type
-      double precision, dimension(:,:), allocatable :: d_compute
-      call lammps_extract_compute_dp2a (d_compute, ptr, id, style, type)
-      allocate (compute(size(d_compute,1), size(d_compute,2)))
-      compute = real(d_compute)
-      deallocate (d_compute)
-   end subroutine lammps_extract_compute_r2a
 
 !-----------------------------------------------------------------------------
 
@@ -616,7 +524,7 @@ contains !! Wrapper functions local to this module {{{1
       fix = lammps_actual_extract_fix (ptr, Cid, Cstyle, Ctype, Ci, Cj)
    end function lammps_extract_fix_Cptr
    subroutine lammps_extract_fix_dp (fix, ptr, id, style, type, i, j)
-      double precision, intent(out) :: fix
+      real (C_double), intent(out) :: fix
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: id
       integer, intent(in) :: style, type, i, j
@@ -635,8 +543,8 @@ contains !! Wrapper functions local to this module {{{1
             call lammps_error_all (ptr, FLERR, 'You cannot extract a fix''s &
                &per-atom/local array (rank 2) into a scalar.')
          case default
-            call lammps_error_all (ptr, FLERR, 'Invalid extract_fix style&
+            call lammps_error_all (ptr, FLERR, 'Invalid extract_fix style/&
-               & value.')
+               &type combination.')
          end select
          return
       end if
@@ -647,22 +555,12 @@ contains !! Wrapper functions local to this module {{{1
       ! Memory is only allocated for "global" fix variables
       if ( style == 0 ) call lammps_free (Cptr)
    end subroutine lammps_extract_fix_dp
-   subroutine lammps_extract_fix_r (fix, ptr, id, style, type, i, j)
-      real, intent(out) :: fix
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: id
-      integer, intent(in) :: style, type, i, j
-      double precision :: d_fix
-      call lammps_extract_fix_dp (d_fix, ptr, id, style, type, i, j)
-      fix = real(d_fix)
-   end subroutine lammps_extract_fix_r
    subroutine lammps_extract_fix_dpa (fix, ptr, id, style, type, i, j)
-      double precision, dimension(:), allocatable, intent(out) :: fix
+      real (C_double), dimension(:), pointer, intent(out) :: fix
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: id
       integer, intent(in) :: style, type, i, j
       type (C_ptr) :: Cptr
-      real (C_double), dimension(:), pointer :: Fptr
       integer :: fix_len
       ! Check for the correct dimensionality
       if ( style == 0 ) then
@@ -682,31 +580,17 @@ contains !! Wrapper functions local to this module {{{1
          return
       end if
       fix_len = lammps_extract_fix_vectorsize (ptr, id, style)
-      allocate (fix(fix_len))
+      call C_F_pointer (Cptr, fix, (/fix_len/))
-      Cptr = lammps_extract_fix_Cptr (ptr, id, style, type, i, j)
+      ! Memory is only allocated for "global" fix variables, which we should
-      call C_F_pointer (Cptr, Fptr, (/fix_len/))
+      ! never get here, so no need to call lammps_free!
-      fix = Fptr
-      nullify (Fptr)
-      ! Memory is only allocated for "global" fix variables
-      if ( style == 0 ) call lammps_free (Cptr)
    end subroutine lammps_extract_fix_dpa
-   subroutine lammps_extract_fix_ra (fix, ptr, id, style, type, i, j)
-      real, dimension(:), allocatable, intent(out) :: fix
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: id
-      integer, intent(in) :: style, type, i, j
-      double precision, dimension(:), allocatable :: d_fix
-      call lammps_extract_fix_dpa (d_fix, ptr, id, style, type, i, j)
-      allocate (fix(size(d_fix)))
-      fix = real(d_fix)
-      deallocate (d_fix)
-   end subroutine lammps_extract_fix_ra
    subroutine lammps_extract_fix_dp2a (fix, ptr, id, style, type, i, j)
-      double precision, dimension(:,:), allocatable, intent(out) :: fix
+      real (C_double), dimension(:,:), pointer, intent(out) :: fix
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: id
       integer, intent(in) :: style, type, i, j
       type (C_ptr) :: Cptr
+      type (C_ptr), pointer, dimension(:) :: Cfix
       integer :: nr, nc
       ! Check for the correct dimensionality
       if ( style == 0 ) then
@@ -723,22 +607,11 @@ contains !! Wrapper functions local to this module {{{1
          return
       end if
       call lammps_extract_fix_arraysize (ptr, id, style, nr, nc)
-      allocate (fix(nr, nc))
+      ! Extract pointer to first element as Cfix(1)
-      Cptr = lammps_extract_fix_Cptr (ptr, id, style, type, i, j)
+      call C_F_pointer (Cptr, Cfix, (/nr/))
-      fix = Cdoublestar_to_2darray (Cptr, nr, nc)
+      ! Now extract the array, which is transposed
-      ! C pointer should not be freed
+      call C_F_pointer (Cfix(1), fix, (/nc, nr/))
    end subroutine lammps_extract_fix_dp2a
-   subroutine lammps_extract_fix_r2a (fix, ptr, id, style, type, i, j)
-      real, dimension(:,:), allocatable, intent(out) :: fix
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: id
-      integer, intent(in) :: style, type, i, j
-      double precision, dimension(:,:), allocatable :: d_fix
-      call lammps_extract_fix_dp2a (d_fix, ptr, id, style, type, i, j)
-      allocate (fix(size(d_fix,1), size(d_fix,2)))
-      fix = real(d_fix)
-      deallocate (d_fix)
-   end subroutine lammps_extract_fix_r2a
 
 !-----------------------------------------------------------------------------
 
@@ -760,24 +633,11 @@ contains !! Wrapper functions local to this module {{{1
       variable = lammps_actual_extract_variable (ptr, Cname, Cgroup)
       deallocate (Cgroup)
    end function lammps_extract_variable_Cptr
-   subroutine lammps_extract_variable_i (variable, ptr, name, group)
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      character (len=*), intent(in), optional :: group
-      integer, intent(out) :: variable
-      double precision :: d_var
-      if ( present (group) ) then
-         call lammps_extract_variable_dp (d_var, ptr, name, group)
-      else
-         call lammps_extract_variable_dp (d_var, ptr, name)
-      end if
-      variable = nint(d_var)
-   end subroutine lammps_extract_variable_i
    subroutine lammps_extract_variable_dp (variable, ptr, name, group)
+      real (C_double), intent(out) :: variable
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: name
       character (len=*), intent(in), optional :: group
-      double precision, intent(out) :: variable
       type (C_ptr) :: Cptr
       real (C_double), pointer :: Fptr
       if ( present(group) ) then
@@ -790,37 +650,8 @@ contains !! Wrapper functions local to this module {{{1
       nullify (Fptr)
       call lammps_free (Cptr)
    end subroutine lammps_extract_variable_dp
-   subroutine lammps_extract_variable_r (variable, ptr, name, group)
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      character (len=*), intent(in), optional :: group
-      real, intent(out) :: variable
-      double precision :: d_var
-      if ( present (group) ) then
-         call lammps_extract_variable_dp (d_var, ptr, name, group)
-      else
-         call lammps_extract_variable_dp (d_var, ptr, name)
-      end if
-      variable = real(d_var)
-   end subroutine lammps_extract_variable_r
-
-   subroutine lammps_extract_variable_ia (variable, ptr, name, group)
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      character (len=*), intent(in), optional :: group
-      integer, dimension(:), allocatable, intent(out) :: variable
-      double precision, dimension(:), allocatable :: d_var
-      if ( present (group) ) then
-         call lammps_extract_variable_dpa (d_var, ptr, name, group)
-      else
-         call lammps_extract_variable_dpa (d_var, ptr, name)
-      end if
-      allocate (variable(size(d_var)))
-      variable = nint(d_var)
-      deallocate (d_var)
-   end subroutine lammps_extract_variable_ia
    subroutine lammps_extract_variable_dpa (variable, ptr, name, group)
-      double precision, dimension(:), allocatable, intent(out) :: variable
+      real (C_double), dimension(:), allocatable, intent(out) :: variable
       type (C_ptr), intent(in) :: ptr
       character (len=*), intent(in) :: name
       character (len=*), intent(in), optional :: group
@@ -839,21 +670,6 @@ contains !! Wrapper functions local to this module {{{1
       nullify (Fptr)
       call lammps_free (Cptr)
    end subroutine lammps_extract_variable_dpa
-   subroutine lammps_extract_variable_ra (variable, ptr, name, group)
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      character (len=*), intent(in), optional :: group
-      real, dimension(:), allocatable, intent(out) :: variable
-      double precision, dimension(:), allocatable :: d_var
-      if ( present (group) ) then
-         call lammps_extract_variable_dpa (d_var, ptr, name, group)
-      else
-         call lammps_extract_variable_dpa (d_var, ptr, name)
-      end if
-      allocate (variable(size(d_var)))
-      variable = real(d_var)
-      deallocate (d_var)
-   end subroutine lammps_extract_variable_ra
 
 !-------------------------------------------------------------------------2}}}
 
@@ -909,17 +725,6 @@ contains !! Wrapper functions local to this module {{{1
       data = Fdata(:)
       deallocate (Fdata)
    end subroutine lammps_gather_atoms_dpa
-   subroutine lammps_gather_atoms_ra (ptr, name, count, data)
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      integer, intent(in) :: count
-      real, dimension(:), allocatable, intent(out) :: data
-      double precision, dimension(:), allocatable :: d_data
-      call lammps_gather_atoms_dpa (ptr, name, count, d_data)
-      allocate (data(size(d_data)))
-      data = d_data
-      deallocate (d_data)
-   end subroutine lammps_gather_atoms_ra
 
 !-----------------------------------------------------------------------------
 
@@ -961,14 +766,6 @@ contains !! Wrapper functions local to this module {{{1
       Cdata = C_loc (Fdata(1))
       call lammps_actual_scatter_atoms (ptr, Cname, Ctype, Ccount, Cdata)
    end subroutine lammps_scatter_atoms_dpa
-   subroutine lammps_scatter_atoms_ra (ptr, name, data)
-      type (C_ptr), intent(in) :: ptr
-      character (len=*), intent(in) :: name
-      real, dimension(:), intent(in) :: data
-      double precision, dimension(size(data)) :: d_data
-      d_data = real (data, kind(d_data))
-      call lammps_scatter_atoms_dpa (ptr, name, d_data)
-   end subroutine lammps_scatter_atoms_ra
 
 !-----------------------------------------------------------------------------
 
@@ -1137,31 +934,6 @@ contains !! Wrapper functions local to this module {{{1
 
    end subroutine Cstring2argcargv
 
-!-----------------------------------------------------------------------------
-
-   function Cdoublestar_to_2darray (Carray, nrows, ncolumns) result (Farray)
-
-   ! Take a C/C++ array of pointers to pointers to doubles (sort of like a
-   ! two-dimensional array, and handled the same way from the programmer's
-   ! perspective) into a Fortran-style array.  Note that columns in C still
-   ! correspond to columns in Fortran here and the same for rows.
-
-      type (C_ptr), intent(in) :: Carray
-      integer, intent(in) :: nrows, ncolumns
-      double precision, dimension(nrows, ncolumns) :: Farray
-      type (C_ptr), dimension(:), pointer :: C_rows
-      real (C_double), dimension(:), pointer :: F_row
-      integer :: i
-
-      ! Convert each "C row pointer" into an array of rows
-      call C_F_pointer (Carray, C_rows, (/nrows/))
-      do i = 1, nrows
-         ! Convert each C pointer (an entire row) into a Fortran pointer
-         call C_F_pointer (C_rows(i), F_row, (/ncolumns/))
-         Farray (i,:) = real(F_row, kind(0.0D0))
-      end do
-
-   end function Cdoublestar_to_2darray
 ! 1}}}
 
 end module LAMMPS

examples/COUPLE/fortran2/README

@@ -4,7 +4,8 @@ src/library.h so they can be used directly from Fortran-encoded programs.
 All functions in src/library.h that use and/or return C-style pointers have
 Fortran wrapper functions that use Fortran-style arrays, pointers, and
 strings; all C-style memory management is handled internally with no user
-intervention.
+intervention.  See --USE-- for notes on how this interface differs from the
+C interface (and the Python interface).
 
 This interface was created by Karl Hammond who you can contact with
 questions:
@@ -25,24 +26,30 @@ You are also advised to read the --USE-- section below before trying to
 compile.
 
 The following steps will work to compile this module (replace ${LAMMPS_SRC}
-with the path to your LAMMPS source directory):
+with the path to your LAMMPS source directory).
- (1) Compile LAMMPS as a static library.  Call the resulting file ${LAMMPS_LIB},
+
-     which will have an actual name lake liblmp_openmpi.a.  If compiling
+Steps 3-5 are accomplished, possibly after some modifications to
-     using the MPI stubs in ${LAMMPS_SRC}/STUBS, you will need to know where
+the makefile, by make using the attached makefile.  Said makefile also builds
-     libmpi.a is as well (I'll call it ${MPI_STUBS} hereafter)
+the dynamically-linkable library (liblammps_fortran.so).
- (2) Copy said library to your Fortran program's source directory or include
+
-     its location in a -L${LAMMPS_SRC} flag to your compiler.
+** STATIC LIBRARY INSTRUCTIONS **
+ (1) Compile LAMMPS as a static library.
+     Call the resulting file ${LAMMPS_LIB}, which will have an actual name
+     like liblmp_openmpi.a.  If compiling using the MPI stubs in
+     ${LAMMPS_SRC}/STUBS, you will need to know where libmpi_stubs.a
+     is as well (I'll call it ${MPI_STUBS} hereafter)
+ (2) Copy said library to your Fortran program's source directory or replace
+     ${LAMMPS_LIB} with its full path in the instructions below.
  (3) Compile (but don't link!) LAMMPS.F90.  Example:
         mpif90 -c LAMMPS.f90
      OR
         gfortran -c LAMMPS.F90
-     Copy the LAMMPS.o and lammps.mod (or whatever your compiler calls module
-     files) to your Fortran program's source directory.
      NOTE:  you may get a warning such as,
         subroutine lammps_open_wrapper (argc, argv, communicator, ptr) &
         Variable 'communicator' at (1) is a parameter to the BIND(C)
         procedure 'lammps_open_wrapper' but may not be C interoperable
      This is normal (see --IMPLEMENTATION NOTES--).
+
  (4) Compile (but don't link) LAMMPS-wrapper.cpp.  You will need its header
      file as well.  You will have to provide the locations of LAMMPS's
      header files.  For example,
@@ -51,13 +58,11 @@ with the path to your LAMMPS source directory):
         g++ -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp
      OR
         icpc -c -I${LAMMPS_SRC} -I${LAMMPS_SRC}/STUBS LAMMPS-wrapper.cpp
-     Copy the resulting object file LAMMPS-wrapper.o to your Fortran program's
+ (5) OPTIONAL:  Make a library from the object files so you can carry around
-     source directory.
+     two files instead of three.  Example:
- (4b) OPTIONAL:  Make a library so you can carry around two files instead of
-     three.  Example:
         ar rs liblammps_fortran.a LAMMPS.o LAMMPS-wrapper.o
      This will create the file liblammps_fortran.a that you can use in place
-     of "LAMMPS.o LAMMPS-wrapper.o" in part (6).  Note that you will still
+     of "LAMMPS.o LAMMPS-wrapper.o" later.  Note that you will still
      need to have the .mod file from part (3).
 
      It is also possible to add LAMMPS.o and LAMMPS-wrapper.o into the
@@ -67,7 +72,7 @@ with the path to your LAMMPS source directory):
      In this case, you can now use the Fortran wrapper functions as if they
      were part of the usual LAMMPS library interface (if you have the module
      file visible to the compiler, that is).
- (5) Compile your Fortran program.  Example:
+ (6) Compile (but don't link) your Fortran program.  Example:
         mpif90 -c myfreeformatfile.f90
         mpif90 -c myfixedformatfile.f
      OR
@@ -78,25 +83,47 @@ with the path to your LAMMPS source directory):
 
      IMPORTANT:  If the Fortran module from part (3) is not in the current
      directory or in one searched by the compiler for module files, you will
-     need to include that location via the -I flag to the compiler.
+     need to include that location via the -I flag to the compiler, like so:
- (6) Link everything together, including any libraries needed by LAMMPS (such
+        mpif90 -I${LAMMPS_SRC}/examples/COUPLE/fortran2 -c myfreeformatfile.f90
+
+ (7) Link everything together, including any libraries needed by LAMMPS (such
      as the C++ standard library, the C math library, the JPEG library, fftw,
      etc.)  For example,
         mpif90 LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
-            ${LAMMPS_LIB} -lstdc++ -lm
+            ${LAMMPS_LIB} -lmpi_cxx -lstdc++ -lm
      OR
         gfortran LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
             ${LAMMPS_LIB} ${MPI_STUBS} -lstdc++ -lm
      OR
         ifort LAMMPS.o LAMMPS-wrapper.o ${my_object_files} \
-            ${LAMMPS_LIB} ${MPI_STUBS} -cxxlib -limf -lm
+            ${LAMMPS_LIB} ${MPI_STUBS} -cxxlib -lm
      Any other required libraries (e.g. -ljpeg, -lfftw) should be added to
      the end of this line.
 
 You should now have a working executable.
 
-Steps 3 and 4 above are accomplished, possibly after some modifications to
+** DYNAMIC LIBRARY INSTRUCTIONS **
-the makefile, by make using the attached makefile.
+ (1) Compile LAMMPS as a dynamic library
+     (make makeshlib && make -f Makefile.shlib [targetname]).
+ (2) Compile, but don't link, LAMMPS.F90 using the -fPIC flag, such as
+        mpif90 -fPIC -c LAMMPS.f90
+ (3) Compile, but don't link, LAMMPS-wrapper.cpp in the same manner, e.g.
+        mpicxx -fPIC -c LAMMPS-wrapper.cpp
+ (4) Make the dynamic library, like so:
+        mpif90 -fPIC -shared -o liblammps_fortran.so LAMMPS.o LAMMPS-wrapper.o
+ (5) Compile your program, such as,
+        mpif90 -I${LAMMPS_SRC}/examples/COUPLE/fortran2 -c myfreeformatfile.f90
+     where ${LAMMPS_SRC}/examples/COUPLE/fortran2 contains the .mod file from
+     step (3)
+ (6) Link everything together, such as
+        mpif90 ${my_object_files} -L${LAMMPS_SRC} \
+            -L${LAMMPS_SRC}/examples/COUPLE/fortran2 -llammps_fortran \
+            -llammps_openmpi -lmpi_cxx -lstdc++ -lm
+
+If you wish to avoid the -L flags, add the directories containing your
+shared libraries to the LIBRARY_PATH environment variable.  At run time, you
+will have to add these directories to LD_LIBRARY_PATH as well; otherwise,
+your executable will not find the libraries it needs.
 
 -------------------------------------
 
@@ -109,8 +136,8 @@ should look something like this:
     ! Other modules, etc.
     implicit none
     type (lammps_instance) :: lmp ! This is a pointer to your LAMMPS instance
-    double precision :: fix
+    real (C_double) :: fix
-    double precision, dimension(:), allocatable :: fix2
+    real (C_double), dimension(:), pointer :: fix2
     ! Rest of declarations
     call lammps_open_no_mpi ('lmp -in /dev/null -screen out.lammps',lmp)
     ! Set up rest of program here
@@ -121,6 +148,11 @@ should look something like this:
   end program call_lammps
 
 Important notes:
+ * Though I dislike the use of pointers, they are necessary when communicating
+   with C and C++, which do not support Fortran's ALLOCATABLE attribute.
+ * There is no need to deallocate C-allocated memory; this is done for you in
+   the cases when it is done (which are all cases when pointers are not
+   accepted, such as global fix data)
  * All arguments which are char* variables in library.cpp are character (len=*)
    variables here.  For example,
         call lammps_command (lmp, 'units metal')
@@ -133,24 +165,27 @@ Important notes:
    as assign a 2D array to a scalar), but it's not perfect.  For example, the
    command
         call lammps_extract_global (nlocal, ptr, 'nlocal')
-   will give nlocal correctly if nlocal is of type INTEGER, but it will give
+   will give nlocal correctly if nlocal is a pointer to type INTEGER, but it
-   the wrong answer if nlocal is of type REAL or DOUBLE PRECISION.  This is a
+   will give the wrong answer if nlocal is a pointer to type REAL.  This is a
    feature of the (void*) type cast in library.cpp.  There is no way I can
-   check this for you!
+   check this for you!  It WILL catch you if you pass it an allocatable or
- * You are allowed to use REAL or DOUBLE PRECISION floating-point numbers.
+   fixed-size array when it expects a pointer.
-   All LAMMPS data (which are of type REAL(C_double)) are rounded off if
+ * Arrays constructed from temporary data from LAMMPS are ALLOCATABLE, and
-   placed in single precision variables.  It is tacitly assumed that NO C++
+   represent COPIES of data, not the originals.  Functions like
-   variables are of type float; everything is int or double (since this is
+   lammps_extract_atom, which return actual LAMMPS data, are pointers.
-   all library.cpp currently handles).
+ * IMPORTANT:  Due to the differences between C and Fortran arrays (C uses
- * An example of a complete program is offered at the end of this file.
+   row-major vectors, Fortran uses column-major vectors), all arrays returned
+   from LAMMPS have their indices swapped.
+ * An example of a complete program, simple.f90, is included with this
+   package.
 
 -------------------------------------
 
 --TROUBLESHOOTING--
 
-Compile-time errors probably indicate that your compiler is not new enough to
+Compile-time errors (when compiling LAMMPS.F90, that is) probably indicate
-support Fortran 2003 features.  For example, GCC 4.1.2 will not compile this
+that your compiler is not new enough to support Fortran 2003 features.  For
-module, but GCC 4.4.0 will.
+example, GCC 4.1.2 will not compile this module, but GCC 4.4.0 will.
 
 If your compiler balks at 'use, intrinsic :: ISO_C_binding,' try removing the
 intrinsic part so it looks like an ordinary module.  However, it is likely
@@ -158,15 +193,15 @@ that such a compiler will also have problems with everything else in the
 file as well.
 
 If you get a segfault as soon as the lammps_open call is made, check that you
-compiled your program AND LAMMPS-header.cpp using the same MPI headers.  Using
+compiled your program AND LAMMPS-wrapper.cpp using the same MPI headers.  Using
-the stubs for one and the actual MPI library for the other will cause major
+the stubs for one and the actual MPI library for the other will cause Bad
-problems.
+Things to happen.
 
 If you find run-time errors, please pass them along via the LAMMPS Users
-mailing list.  Please provide a minimal working example along with the names
+mailing list (please CC me as well; address above).  Please provide a minimal
-and versions of the compilers you are using.  Please make sure the error is
+working example along with the names and versions of the compilers you are
-repeatable and is in MY code, not yours (generating a minimal working example
+using.  Please make sure the error is repeatable and is in MY code, not yours
-will usually ensure this anyway).
+(generating a minimal working example will usually ensure this anyway).
 
 -------------------------------------
 
@@ -177,22 +212,23 @@ their purpose is the same, but they may take different arguments.  Here are
 some of the important differences:
  * lammps_open and lammps_open_no_mpi take a string instead of argc and
    argv.  This is necessary because C and C++ have a very different way
-   of treating strings than Fortran.
+   of treating strings than Fortran.  If you want the command line to be
+   passed to lammps_open (as it often would be from C/C++), use the
+   GET_COMMAND intrinsic to obtain it.
  * All C++ functions that accept char* pointers now accept Fortran-style
    strings within this interface instead.
  * All of the lammps_extract_[something] functions, which return void*
    C-style pointers, have been replaced by generic subroutines that return
    Fortran variables (which may be arrays).  The first argument houses the
-   variable to be returned; all other arguments are identical except as
+   variable/pointer to be returned (pretend it's on the left-hand side); all
-   stipulated above.  Note that it is not possible to declare generic
+   other arguments are identical except as stipulated above.
-   functions that are selected based solely on the type/kind/rank (TKR)
+   Note that it is not possible to declare generic functions that are selected
-   signature of the return value, only based on the TKR of the arguments.
+   based solely on the type/kind/rank (TKR) signature of the return value,
+   only based on the TKR of the arguments.
  * The SHAPE of the first argument to lammps_extract_[something] is checked
    against the "shape" of the C array (e.g., double vs. double* vs. double**).
    Calling a subroutine with arguments of inappropriate rank will result in an
    error at run time.
- * All arrays passed to subroutines must be ALLOCATABLE and are REALLOCATED
-   to fit the shape of the array LAMMPS will be returning.
  * The indices i and j in lammps_extract_fix are used the same way they
    are in f_ID[i][j] references in LAMMPS (i.e., starting from 1).  This is
    different than the way library.cpp uses these numbers, but is more
@@ -202,8 +238,7 @@ some of the important differences:
    instead of a function.
  * The pointer to LAMMPS itself is of type(lammps_instance), which is itself
    a synonym for type(C_ptr), part of ISO_C_BINDING.  Type (C_ptr) is
-   C's void* data type.  This should be the only C data type that needs to
+   C's void* data type.
-   be used by the end user.
  * This module will almost certainly generate a compile-time warning,
    such as,
          subroutine lammps_open_wrapper (argc, argv, communicator, ptr) &
@@ -213,9 +248,18 @@ some of the important differences:
    INTEGER argument, whose type is defined by the MPI library itself.  The
    Fortran integer is converted to a C integer by the MPI library (if such
    conversion is actually necessary).
- * Unlike library.cpp, this module returns COPIES of the data LAMMPS actually
+  * lammps_extract_global returns COPIES of the (scalar) data, as does the
-   uses.  This is done for safety reasons, as you should, in general, not be
+    C version.
-   overwriting LAMMPS data directly from Fortran.  If you require this
+  * lammps_extract_atom, lammps_extract_compute, and lammps_extract_fix
-   functionality, it is possible to write another function that, for example,
+        have a first argument that will be associated with ACTUAL LAMMPS DATA.
-   returns a Fortran pointer that resolves to the C/C++ data instead of
+        This means the first argument must be:
-   copying the contents of that pointer to the original array as is done now.
+          * The right rank (via the DIMENSION modifier)
+          * A C-interoperable POINTER type (i.e., INTEGER (C_int) or
+            REAL (C_double)).
+  * lammps_extract_variable returns COPIES of the data, as the C library
+    interface does.  There is no need to deallocate using lammps_free.
+  * The 'data' argument to lammps_gather_atoms and lammps_scatter atoms must
+    be ALLOCATABLE.  It should be of type INTEGER or DOUBLE PRECISION.  It
+    does NOT need to be C inter-operable (and indeed should not be).
+  * The 'count' argument of lammps_scatter_atoms is unnecessary; the shape of
+    the array determines the number of elements LAMMPS will read.

examples/COUPLE/fortran2/in.simple

@@ -1,10 +1,11 @@
-units metal
+units lj
-lattice bcc 3.1656
+atom_modify map array
+lattice bcc 1.0
 region simbox block 0 10 0 10 0 10
 create_box 2 simbox
 create_atoms 1 region simbox
-pair_style eam/fs
+pair_style lj/cut 2.5
-pair_coeff * * path/to/my_potential.eam.fs A1 A2
+pair_coeff * * 1.0 1.0
 mass 1 58.2 # These are made-up numbers
 mass 2 28.3
 velocity all create 1200.0 7474848 dist gaussian

examples/COUPLE/fortran2/simple.f90

@@ -1,44 +1,111 @@
 program simple
 
+   use MPI
    use LAMMPS
+
+   ! The following line is unnecessary, as I have included these three entities
+   ! with the LAMMPS module, but I leave them in anyway to remind people where
+   ! they came from
+   use, intrinsic :: ISO_C_binding, only : C_double, C_ptr, C_int
+
    implicit none
 
-   type (lammps_instance) :: lmp
+   ! Notes:
-   double precision :: compute, fix, fix2
+   !  * If LAMMPS returns a scalar that is allocated by the library interface
-   double precision, dimension(:), allocatable :: compute_v, mass, r
+   !     (see library.cpp), then that memory is deallocated automatically and
-   double precision, dimension(:,:), allocatable :: x
+   !     the argument to lammps_extract_fix must be a SCALAR.
-   real, dimension(:,:), allocatable :: x_r
+   !  * If LAMMPS returns a pointer to an array, consisting of internal LAMMPS
+   !     data, then the argument must be an interoperable Fortran pointer.
+   !     Interoperable means it is of type INTEGER (C_INT) or of type
+   !     REAL (C_DOUBLE) in this context.
+   !  * Pointers should NEVER be deallocated, as that would deallocate internal
+   !     LAMMPS data!
+   !  * Note that just because you can read the values of, say, a compute at
+   !     any time does not mean those values represent the "correct" values.
+   !     LAMMPS will abort you if you try to grab a pointer to a non-current
+   !     entity, but once it's bound, it's your responsibility to check that
+   !     it's current before evaluating.
+   !  * IMPORTANT:  Two-dimensional arrays (such as 'x' from extract_atom)
+   !     will be transposed from what they might look like in C++.  This is
+   !     because of different bookkeeping conventions between Fortran and C
+   !     that date back to about 1970 or so (when C was written).
+   !  * Arrays start from 1, EXCEPT for mass from extract_atom, which
+   !     starts from 0.  This is because the C array actually has a blank
+   !     first element (and thus mass[1] corresponds to the mass of type 1)
+
+   type (C_ptr) :: lmp
+   real (C_double), pointer :: compute => NULL()
+   real (C_double) :: fix, fix2
+   real (C_double), dimension(:), pointer :: compute_v => NULL()
+   real (C_double), dimension(:,:), pointer :: x => NULL()
+   real (C_double), dimension(:), pointer :: mass => NULL()
+   integer, dimension(:), allocatable :: types
+   double precision, dimension(:), allocatable :: r
+   integer :: error, narg, me, nprocs
+   character (len=1024) :: command_line
+
+   call MPI_Init (error)
+   call MPI_Comm_rank (MPI_COMM_WORLD, me, error)
+   call MPI_Comm_size (MPI_COMM_WORLD, nprocs, error)
+
+   ! You are free to pass any string you like to lammps_open or
+   ! lammps_open_no_mpi; here is how you pass it the command line
+   !call get_command (command_line)
+   !call lammps_open (command_line, MPI_COMM_WORLD, lmp)
+
+   ! And here's how to to it with a string constant of your choice
+   call lammps_open_no_mpi ('lmp -log log.simple', lmp)
 
-   call lammps_open_no_mpi ('',lmp)
    call lammps_file (lmp, 'in.simple')
    call lammps_command (lmp, 'run 500')
 
+   ! This extracts f_2 as a scalar (the last two arguments can be arbitrary)
    call lammps_extract_fix (fix, lmp, '2', 0, 1, 1, 1)
    print *, 'Fix is ', fix
 
+   ! This extracts f_4[1][1] as a scalar
    call lammps_extract_fix (fix2, lmp, '4', 0, 2, 1, 1)
    print *, 'Fix 2 is ', fix2
 
+   ! This extracts the scalar compute of compute thermo_temp
    call lammps_extract_compute (compute, lmp, 'thermo_temp', 0, 0)
    print *, 'Compute is ', compute
 
+   ! This extracts the vector compute of compute thermo_temp
    call lammps_extract_compute (compute_v, lmp, 'thermo_temp', 0, 1)
    print *, 'Vector is ', compute_v
 
+   ! This extracts the masses
    call lammps_extract_atom (mass, lmp, 'mass')
-   print *, 'Mass is ', mass
+   print *, 'Mass is ', mass(1:)
 
+   ! Extracts a pointer to the arrays of positions for all atoms
    call lammps_extract_atom (x, lmp, 'x')
-   if ( .not. allocated (x) ) print *, 'x is not allocated'
+   if ( .not. associated (x) ) print *, 'x is not associated'
-   print *, 'x is ', x(1,:)
+   print *, 'x is ', x(:,1)  ! Prints x, y, z for atom 1
 
-   call lammps_extract_atom (x_r, lmp, 'x')
+   ! Extracts pointer to atom types
-   if ( .not. allocated (x_r) ) print *, 'x is not allocated'
+   call lammps_gather_atoms (lmp, 'type', 1, types)
-   print *, 'x_r is ', x_r(1,:)
+   print *, 'types is ', types(1:3)
 
-   call lammps_get_coords (lmp, r)
+   ! Allocates an array and assigns all positions to it
-   print *, 'r is ', r(1:3)
+   call lammps_gather_atoms (lmp, 'x', 3, r)
+   print *, 'size(r) = ', size(r)
+   print *, 'r is ', r(1:6)
+
+   ! Puts those position data back
+   call lammps_scatter_atoms (lmp, 'x', r)
+
+   call lammps_command (lmp, 'run 1')
+   print *, 'x is ', x(:,1)  ! Note that the position updates!
+   print *, 'Compute is ', compute  ! This did only because "temp" is part of
+                                    ! the thermo output; the vector part did
+                                    ! not, and won't until we give LAMMPS a
+                                    ! thermo output or other command that
+                                    ! requires its value
 
    call lammps_close (lmp)
 
+   call MPI_Finalize (error)
+
 end program simple