compute_grid Subroutine

subroutine compute_grid(mol, cutoff, grid, det, lat_inv, cell_w)

   !> Stucture type
   type(structure_type), intent(in) :: mol

   !> Interaction cutoff radius
   real(wp), intent(in) :: cutoff

   !> Linked cell grid type
   type(grid_type), intent(out) :: grid

   !> Determinant (Volume) of the lattice
   real(wp), intent(out) :: det

   !> Inverse of the lattice matrix, zero for a molecular system
   real(wp), intent(out), optional :: lat_inv(3, 3)

   !> Perpendicular width of each grid cell
   real(wp), intent(out), optional :: cell_w(3)

   real(wp) :: H(3), cross_ij(3), lattice(3, 3)
   real(wp) :: max_xyz(3), min_xyz(3)
   integer  :: i

   if (present(lat_inv)) lat_inv = 0.0_wp

   if (any(mol%periodic)) then
      ! Inverse Lattice Matrix
      lattice = mol%lattice
      det = matdet_3x3(lattice)

      if (present(lat_inv)) lat_inv = matinv_3x3(lattice)

      ! Perpendicular height H1 (normal to a2 x a3)
      cross_ij(:) = crossprod(lattice(:,2), lattice(:,3))
      H(1) = abs(det) / sqrt(sum(cross_ij**2))

      ! Perpendicular height H2 (normal to a3 x a1)
      cross_ij(:) = crossprod(lattice(:,3), lattice(:,1))
      H(2) = abs(det) / sqrt(sum(cross_ij**2))

      ! Perpendicular height H3 (normal to a1 x a2)
      cross_ij(:) = crossprod(lattice(:,1), lattice(:,2))
      H(3) = abs(det) / sqrt(sum(cross_ij**2))

      ! Number of cells: must be at least 1, and linked cell height >= cutoff
      do i = 1, 3
         grid%n_xyz(i) = max(1, floor(H(i) / (cutoff + tiny(1.0_wp))))
      end do

      if (present(cell_w)) cell_w = H / real(grid%n_xyz, wp)
   else
      min_xyz = minval(mol%xyz, dim=2) - buffer
      max_xyz = maxval(mol%xyz, dim=2) + buffer

      ! Number of cells: must be at least 1, and cell width >= cutoff
      grid%n_xyz = max(1, floor((max_xyz - min_xyz) / (cutoff + tiny(1.0_wp))))

      if (present(cell_w)) then
         cell_w = (max_xyz - min_xyz) / (real(grid%n_xyz, wp) + tiny(1.0_wp)) &
            & + tiny(1.0_wp)
      end if
      det = product(max_xyz - min_xyz)
   end if

   grid%ncell = grid%n_xyz(1)*grid%n_xyz(2)*grid%n_xyz(3)
   grid%periodic = any(mol%periodic)

end subroutine compute_grid