analyze_wavefunction_symmetry Subroutine

    subroutine analyze_wavefunction_symmetry()
        ! driver routine to analyze the symmetry of a part of the
        ! wavefunction.
        integer :: sym_orbs(nBasis / 2), orig_orbs(nBasis / 2)
        integer :: i, transformed_states(nel, n_symmetry_states)
        real(dp) :: transformed_weights(n_symmetry_states)
        integer(n_int) :: transformed_states_ilut(0:niftot, n_symmetry_states)

        ! i have the original list of the orbitals and the
        ! corresponding k- or r-vectors.

        call init_symmetry_states()

        if (iProcIndex == root) then
            print *, "Analyze the symmetry of the FCIQMC wavefunction: "

            if (t_symmetry_rotation) then
                print *, "applying rotation of ", symmetry_rotation_angle, " degrees"
            end if

            if (t_symmetry_mirror) then
                print *, "mirror wf. along ", symmertry_mirror_axis, "-axis"
            end if

            if (t_symmetry_inversion) then
                print *, "apply inversion symmetry"
            end if

            if (t_read_symmetry_states) then
                print *, "using the ", n_symmetry_states, " read-in states: "
                do i = 1, n_symmetry_states
                    print *, symmetry_states(:, i)
                end do
            end if

            if (t_pop_symmetry_states) then
                print *, "using the ", n_symmetry_states, " highest occupied states"
                do i = 1, n_symmetry_states
                    print *, symmetry_states(:, i)
                end do
            end if

            ! first check if everything is setup correctly..
            print *, "brr:"
            do i = 1, nBasis
                print *, brr(i)
            end do

            print *, "i, lat%ind, lat%get_sym, lat%get_k_vec"
            do i = 1, nBasis / 2
                print *, i, lat%get_site_index(i), lat%get_sym(i), lat%get_k_vec(i)
            end do

            call writebasis(stdout, g1, nBasis, arr, brr)

            ! i need to apply the chosen symmetry to the vectors and determine,
            ! which orbital maps into which
            ! the flag if we actually apply is within the functions.
            ! ROTATION:
            orig_orbs = get_spatial(brr(1:nBasis:2))
            sym_orbs = apply_rotation(orig_orbs, symmetry_rotation_angle)

            print *, "orig orbs -> sym_orbs: "
            do i = 1, nBasis / 2
                print *, orig_orbs(i), " -> ", sym_orbs(i)
            end do

            ! i can ofc apply multiple symmetries
            ! MIRROR
            sym_orbs = apply_mirror(sym_orbs, symmertry_mirror_axis)

            ! INVERSION
            sym_orbs = apply_inversion(sym_orbs)

            ! now i have the mapping between the original and the
            ! symmetry transformed orbitals

            call transform_states(orig_orbs, sym_orbs, n_symmetry_states, &
                                  symmetry_states, symmetry_weights, symmetry_states_ilut, transformed_states, &
                                  transformed_weights, transformed_states_ilut)

            call compare_states(n_symmetry_states, symmetry_states_ilut, transformed_states_ilut)

        end if

    end subroutine analyze_wavefunction_symmetry