VASPBasisInit Subroutine

subroutine VASPBasisInit(ARR, BRR, G1, LEN)
    ! Largely lifted from the CPMD analogue.  Should be doing (roughly) the same
    ! thing to get going!
    use SystemData, only: Symmetry, SymmetrySize, SymmetrySizeB
    use SystemData, only: BasisFN, BasisFNSize, BasisFNSizeB, nBASISMax, NullBasisFn
    use vasp_interface, only: nStates, nKP, KPntInd, eigv
    use SymData, only: KPntSym, nSym
    use constants, only: dp, stdout
    use sym_mod, only: ComposeAbelianSym, GenKPtIrreps
    implicit none
    integer :: LEN
    real(dp) :: ARR(LEN, 2)
    integer :: BRR(LEN)
    type(BasisFN) :: G1(LEN)
    integer :: i
    type(Symmetry) :: iDecomp

    NBASISMAX(1:3, 1:2) = 0
    NBASISMAX(1, 3) = 2
    NBASISMAX(4, 1) = -1
    NBASISMAX(4, 2) = 1

!  set ISPINSKIP=0, to tell the SCRs that there's no UMAT
    NBASISMAX(2, 3) = 0

    call GenKPtIrreps(NKP, NKP, KPNTIND, NSTATES)

    G1(1:LEN) = NullBasisFn

    do i = 1, nStates
        IDECOMP%s = ComposeAbelianSym(KpntSym(:, KPntInd(I)))
        G1(I * 2 - 1)%SYM = IDECOMP
        G1(I * 2)%SYM = IDECOMP
        G1(I * 2 - 1)%MS = -1
        G1(I * 2)%MS = 1
        G1(I * 2 - 1)%K(1:3) = 0
        G1(I * 2)%K(1:3) = 0
        ARR(2 * I - 1, 1) = eigv(I)
        ARR(2 * I, 1) = eigv(I)
        ARR(2 * I - 1, 2) = eigv(I)
        ARR(2 * I, 2) = eigv(I)
        BRR(2 * I - 1) = 2 * I - 1
        BRR(2 * I) = 2 * I
    end do

    write(stdout, *) 'Using Abelian symmetry formulation.'
    NBASISMAX(5, 1) = 0
    NBASISMAX(5, 2) = NSYM - 1

! Show it's a generic spatial basis
    NBASISMAX(3, 3) = 1

    return
end subroutine VASPBasisInit