diff --git a/src/crystalmod@env.f90 b/src/crystalmod@env.f90
index 6eabdca5..0df03a50 100644
--- a/src/crystalmod@env.f90
+++ b/src/crystalmod@env.f90
@@ -74,7 +74,8 @@ end subroutine build_env
   !> - up2dcidx = up to a distance to each atom in the complete list.
   !> - up2sh = up to a number of shells.
   !> - up2n = up to a number of atoms.
-  !> The atom list with up2n is sorted by distance on output.
+  !> If up2n, up2sh, or ishell0 are present, the atom list is sorted by distance on output
+  !> regardless of the value of sorted.
   !>
   !> Optional input:
   !> - nid0 = consider only atoms with index nid0 from the nneq list.
@@ -139,7 +140,7 @@ module subroutine list_near_atoms(c,xp,icrd,sorted,nat,eid,dist,lvec,ishell0,up2
 
     ! process input
     nozero_ = .false.
-    sorted_ = sorted .or. present(up2n) .or. present(up2sh)
+    sorted_ = sorted .or. present(up2n) .or. present(up2sh) .or. present(ishell0)
     if (present(nozero)) nozero_ = nozero
     if (present(up2d)) then
        dmax = up2d
@@ -318,7 +319,7 @@ module subroutine list_near_atoms(c,xp,icrd,sorted,nat,eid,dist,lvec,ishell0,up2
              else
                 dmax = rshel(iord(up2sh)) + shell_eps
              end if
-             deallocate(rshel,iord)
+             deallocate(iord)
           end if
 
           ! filter out the unneeded atoms
@@ -438,6 +439,21 @@ module subroutine list_near_atoms(c,xp,icrd,sorted,nat,eid,dist,lvec,ishell0,up2
        deallocate(iord)
     end if
 
+    ! assign shells if requested
+    if (present(ishell0)) then
+       if (allocated(ishell0)) deallocate(ishell0)
+       allocate(ishell0(nat))
+       dd = -1d0
+       k = 0
+       do i = 1, nat
+          if (abs(at_dist(i) - dd) > shell_eps) then
+             k = k + 1
+             dd = at_dist(i)
+          end if
+          ishell0(i) = k
+       end do
+    end if
+
     ! reduce the list if up2n (always sorted)
     if (present(up2n)) nat = up2n_
 
diff --git a/src/global@proc.F90 b/src/global@proc.F90
index fefc110a..f2e8ef09 100644
--- a/src/global@proc.F90
+++ b/src/global@proc.F90
@@ -885,8 +885,8 @@ module subroutine critic_setvariables(line,lp)
     real*8 :: x(3), xx(3), dmax ! xxxx
     integer :: i, j, nat, ierr, nat1, nat2 ! xxxx
     integer :: nid, lvec(3) ! xxxx
-    integer, allocatable :: eid(:), lvecx(:,:)
-    real*8, allocatable :: dist(:)
+    integer, allocatable :: eid(:), lvecx(:,:), ishell0(:) ! xxxx
+    real*8, allocatable :: dist(:) ! xxxx
     real*8 :: dist ! xxxx
     real*8 :: f ! xxxx
     real*8 :: fp(3) ! xxxx
@@ -1086,11 +1086,11 @@ module subroutine critic_setvariables(line,lp)
        ! x = (/0.9d0,0.1d0,0.3d0/)
        ! x=(/0.151912d0,0.20933d0,0.00913d0/)
 
-       call tictac("")
-       do i = 1, 1000
-          call sy%c%env%list_near_atoms(x,icrd_crys,.true.,nat,ierr,eid,dist,lvec,up2n=10)
-       end do
-       call tictac("")
+       ! call tictac("")
+       ! do i = 1, 1000
+       !    call sy%c%env%list_near_atoms(x,icrd_crys,.true.,nat,ierr,eid,dist,lvec,up2n=10)
+       ! end do
+       ! call tictac("")
        ! call sy%c%env%list_near_atoms(x,icrd_crys,.true.,nat,ierr,eid,dist,lvec,up2n=10)
        ! ! call sy%c%env%promolecular(x,icrd_crys,f,fp,fpp,2)
        ! do i = 1, nat
@@ -1103,21 +1103,22 @@ module subroutine critic_setvariables(line,lp)
        ! end do
        ! write (*,*)
 
-       call tictac("")
-       do i = 1, 1000
-          call sy%c%list_near_atoms(x,icrd_crys,.true.,nat,eid,dist,lvec=lvecx,up2n=10)
-       end do
-       call tictac("")
-       ! call sy%c%list_near_atoms(x,icrd_crys,.true.,nat,eid,dist,lvec=lvecx,up2sh=10)
-       ! ! call sy%c%promolecular_env(x,icrd_crys,f,fp,fpp,2)
-       ! do i = 1, nat
-       !    xx = sy%c%atcel(eid(i))%x + lvecx(:,i)
-       !    if (sy%c%ismolecule) then
-       !       xx = (sy%c%x2c(xx) + sy%c%molx0) * 0.529177d0
-       !    end if
-       !    write (*,*) i, xx, dist(i)
+       ! call tictac("")
+       ! do i = 1, 1000
+       !    call sy%c%list_near_atoms(x,icrd_crys,.true.,nat,eid,dist,lvec=lvecx,up2n=10)
        ! end do
-       ! write (*,*)
+       ! call tictac("")
+
+       call sy%c%list_near_atoms(x,icrd_crys,.true.,nat,eid,dist,lvec=lvecx,&
+          ishell0=ishell0,up2d=10d0)
+       do i = 1, nat
+          xx = sy%c%atcel(eid(i))%x + lvecx(:,i)
+          if (sy%c%ismolecule) then
+             xx = (sy%c%x2c(xx) + sy%c%molx0) * 0.529177d0
+          end if
+          write (*,*) i, xx, dist(i), ishell0(i)
+       end do
+       write (*,*)
 
        stop 1