setup_momentum_operators Subroutine

    subroutine setup_momentum_operators()
        use SystemData, only: G1
        implicit none

        integer :: k(3), state, momentum_state, spin, i, j

        if (gf_type == -1) then
            momentum_state = pops_pert(1)%ann_orbs(1)
        else if (gf_type == 1) then
            momentum_state = pops_pert(1)%crtn_orbs(1)
        else
            call stop_all("Equalize initial phase", "Invalid momentum space green's function")
        end if
        allocate(phase_factors(nBasis / 2))
        k = G1(momentum_state)%k
        spin = G1(momentum_state)%ms
        call clear_pops_pert(pops_pert)
        call clear_pops_pert(overlap_pert)
        allocate(pops_pert(nBasis / 2))
        allocate(overlap_pert(nBasis / 2))
        i = 0
        do state = 1, nBasis
            if (G1(state)%ms /= spin) cycle
            i = i + 1
            phase_factors(i) = 0
            do j = 1, 3
                phase_factors(i) = phase_factors(i) + G1(state)%k(j) * k(j)
            end do
            call setup_single_perturbation(i, state, gf_type)
        end do
    end subroutine setup_momentum_operators