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:                             ! Have found an allowed wavevector.
                              ! Add 2 spin orbitals to the set of the basis functions.
                              ibasis = ibasis + 1
                              call init_kpoint(basis_fns(ibasis), kp(1:ndim), 1)
                              ibasis = ibasis + 1
                              call init_kpoint(basis_fns(ibasis), kp(1:ndim), -1)
                          end if
                      end if
                  end do
  
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
  
  