print_d4h_pg Subroutine

    subroutine print_d4h_pg(states)
        ! construct the d4h symmetry operation matrix representations
        ! the operation for now are: E, 2C4, C2, Mv, Md, and inversion
        integer, intent(in) :: states(:, :)

        integer :: i, orig_orbs(nBasis / 2), trans_orbs(nBasis / 2), &
                   matrix_rep(size(states, 2), size(states, 2)), &
                   temp_states(nel, size(states, 2)), phase

        orig_orbs = get_spatial(brr(1:nBasis:2))

        ! first E:
        trans_orbs = orig_orbs
        matrix_rep = construct_matrix_representation(states, orig_orbs, trans_orbs)
        print *, "E:"
        call print_matrix(matrix_rep)
        print *, "character: ", sum([(matrix_rep(i, i), i=1, size(matrix_rep, 1))])

        ! first C4:
        trans_orbs = apply_rotation(orig_orbs, 90.0_dp)
        matrix_rep = construct_matrix_representation(states, orig_orbs, trans_orbs)
        print *, "C4: "
        call print_matrix(matrix_rep)
        print *, "character: ", sum([(matrix_rep(i, i), i=1, size(matrix_rep, 1))])

!         ! first C4:
!         trans_orbs = apply_rotation(orig_orbs, 270.0)
!         matrix_rep = construct_matrix_representation(states, orig_orbs, trans_orbs)
!         print *, "2C4: "
!         call print_matrix(matrix_rep)
!         print *, "character: ", sum([(matrix_rep(i,i), i = 1, size(matrix_rep,1))])

        ! C2:
        trans_orbs = apply_rotation(orig_orbs, 180.0_dp)
        matrix_rep = construct_matrix_representation(states, orig_orbs, trans_orbs)
        print *, "C2: "
        call print_matrix(matrix_rep)
        print *, "character: ", sum([(matrix_rep(i, i), i=1, size(matrix_rep, 1))])

        ! Mv
        trans_orbs = apply_mirror(orig_orbs, 'x')
        matrix_rep = construct_matrix_representation(states, orig_orbs, trans_orbs)
        print *, "Mx: "
        call print_matrix(matrix_rep)
        print *, "character: ", sum([(matrix_rep(i, i), i=1, size(matrix_rep, 1))])

        ! Md
        trans_orbs = apply_mirror(orig_orbs, 'd')
        matrix_rep = construct_matrix_representation(states, orig_orbs, trans_orbs)
        print *, "Md: "
        call print_matrix(matrix_rep)
        print *, "character: ", sum([(matrix_rep(i, i), i=1, size(matrix_rep, 1))])

        ! i
        trans_orbs = apply_inversion(orig_orbs)
        matrix_rep = construct_matrix_representation(states, orig_orbs, trans_orbs)
        print *, "i: "
        call print_matrix(matrix_rep)
        print *, "character: ", sum([(matrix_rep(i, i), i=1, size(matrix_rep, 1))])

    end subroutine print_d4h_pg