anchor.diff:          nmax = 0 ! Set nmax(i) to be 0 for unused higher dimensions.
           forall (i=1:ndim) nmax(i) = maxval(abs(nint(box_length(i)**2/(2*lattice(:,i)))))
   
  -        allocate(basis_fns(2*nsites), stat=ierr)
  +        allocate(basis_fns(nbasis), stat=ierr)
  +        allocate(tmp_basis_fns(nbasis), stat=ierr)
  +        allocate(basis_fns_ranking(nbasis), stat=ierr)
   
           ibasis = 0
           ! The following can't be done as a forall due to compilers being crap...
  
positive.code: 
          if (ibasis /= nbasis) call stop_all('init_basis_fns','Not enough basis functions found.')
  
          if (iproc == parent) then
              write (6,'(/,1X,a15,/,1X,15("-"),/)') 'Basis functions'
              write (6,'(1X,a7)', advance='no') 'k-point'
              do i = 1,ndim
                  write (6,'(3X)', advance='no')
              end do
              write (6,'(a,4X,a7)') 'ms','kinetic'
              do i = 1,nbasis
                  call write_kpoint(basis_fns(i))
              end do
          end if
  
negative.code:         ! the list of wavevectors and hence the kinetic energy associated
          ! with each basis function (two per wavevector to account for spin).
  
          use m_mrgref, only: mrgref
          use errors, only: stop_all
          use parallel, only: iproc, parent
  
          integer :: nmax(3), kp(3) ! Support a maximum of 3 dimensions.
          integer :: i, j, k, ibasis, ierr
          type(kpoint), allocatable, target :: tmp_basis_fns(:)
          type(kpoint), pointer :: basis_fn
          integer, allocatable :: basis_fns_ranking(:)
  
          nbasis = 2*nsites
  
  