get_lattice_points Subroutine

subroutine get_lattice_points(periodic, lat, rthr, trans)
   logical, intent(in) :: periodic(:)
   real(wp), intent(in) :: rthr
   real(wp), intent(in) :: lat(:, :)
   real(wp), allocatable, intent(out) :: trans(:, :)
   real(wp) :: vec(size(lat, 2)), norms(size(lat, 1), size(lat, 2))
   integer :: rep(3)
   integer :: itr, ix, iy, iz

   if (.not.any(periodic)) then
      allocate(trans(3, 1))
      trans(:, :) = 0.0_wp
      return
   end if

   if (all(periodic)) then
      call get_translations(lat, rthr, rep)
   else if (count(periodic) == 1) then
      vec = sqrt(sum(lat**2, 1))
      where(periodic)
         rep = ceiling(rthr / vec)
      elsewhere
         rep = 0
      endwhere
   else if (count(periodic) == 2) then
      call get_normals(lat, norms)
      where(spread(periodic, 2, 3))
         norms = lat
      endwhere
      call get_translations(norms, rthr, rep)
      where(.not.periodic)
         rep = 0
      endwhere
   end if

   allocate(trans(3, product(2*rep+1)))
   itr = 0
   do ix = -rep(1), rep(1)
      do iy = -rep(2), rep(2)
         do iz = -rep(3), rep(3)
            itr = itr + 1
            trans(:, itr) = lat(:, 1)*ix + lat(:, 2)*iy + lat(:, 3)*iz
         end do
      end do
   end do

end subroutine get_lattice_points