setup_kPointToBasisFn Subroutine

    subroutine setup_kPointToBasisFn(in_lat)
        class(lattice), intent(in), optional :: in_lat
        character(*), parameter :: this_routine = "setup_kPointToBasisFn"

        integer :: i, kmaxX, kminX, kmaxY, kminY, kminZ, kmaxZ, iSpinIndex

        if (present(in_lat)) then

            if (allocated(kPointToBasisFn)) return

            if (all(nBasisMax == 0)) then
                call setup_nbasismax(in_lat)
                call setup_g1(in_lat)
            end if

            if (.not. associated(G1)) then
                call setup_g1(in_lat)
            end if

            kmaxX = 0
            kminX = 0
            kmaxY = 0
            kminY = 0
            ! [W.D:] can we make this the same as in the UEG:
            kminZ = 0
            kmaxZ = 0
            do i = 1, 2 * in_lat%get_nsites()
                ! In the hubbard model with tilted lattice boundary conditions,
                ! it's unobvious what the maximum values of
                ! kx and ky are, so this should be found
                IF (G1(i)%k(1) > kmaxX) kmaxX = G1(i)%k(1)
                IF (G1(i)%k(1) < kminX) kminX = G1(i)%k(1)
                IF (G1(i)%k(2) > kmaxY) kmaxY = G1(i)%k(2)
                IF (G1(i)%k(2) < kminY) kminY = G1(i)%k(2)
                if (G1(i)%k(3) > kmaxz) kmaxz = g1(i)%k(3)
                if (G1(i)%k(3) < kminz) kminz = g1(i)%k(3)
            end do
            allocate(kPointToBasisFn(kminX:kmaxX, kminY:kmaxY, kminz:kmaxz, 2))
            !Init to invalid
            kPointToBasisFn = -1
            do i = 1, 2 * in_lat%get_nsites()
                ! iSpinIndex equals 1 for a beta spin (ms=-1), and 2 for an alpha spin (ms=1)
                iSpinIndex = (G1(i)%Ms + 1) / 2 + 1
                kPointToBasisFn(G1(i)%k(1), G1(i)%k(2), G1(i)%k(3), iSpinIndex) = i
            end do
        else
            call Stop_All(this_routine, "not yet implemented!")
        end if

    end subroutine setup_kPointToBasisFn