setup_ras_class Subroutine

    subroutine setup_ras_class(ras, ras_class)

        ! Create a one-to-one map from the original address of a string, as found by
        ! the function get_address, to a new address where states are sorted in order of
        ! their symmetry. Also count the number of strings with each symmetry label.

        type(ras_parameters), intent(in) :: ras
        type(ras_class_data), intent(inout) :: ras_class
        integer :: string(ras_class%nelec_1 + ras_class%nelec_2 + ras_class%nelec_3)
        integer, allocatable, dimension(:) :: symmetries, addresses
        integer :: i, j, counter
        logical :: none_left

        ! Find the number of strings in this class, the maximum address. To find this, put
        ! all electrons in their maximum possible orbtials.
        counter = 0
        do j = 1, ras_class%nelec_1
            counter = counter + 1
            string(counter) = ras%size_1 - ras_class%nelec_1 + j
        end do
        do j = 1, ras_class%nelec_2
            counter = counter + 1
            string(counter) = ras%size_1 + ras%size_2 - ras_class%nelec_2 + j
        end do
        do j = 1, ras_class%nelec_3
            counter = counter + 1
            string(counter) = ras%size_1 + ras%size_2 + ras%size_3 - ras_class%nelec_3 + j
        end do

        ras_class%class_size = get_address(ras_class, string)

        ! Now that we have the class size.
        allocate(symmetries(ras_class%class_size))
        allocate(ras_class%address_map(ras_class%class_size))
        allocate(addresses(ras_class%class_size))

        ! Now generate the string with the lowest address.
        call generate_first_full_string(string, ras, ras_class)
        symmetries(1) = get_abelian_sym(string)
        addresses(1) = 1

        ! Loop over all possible strings in this class, and store the symmetry of each.
        do i = 2, ras_class%class_size
            call generate_next_string(string, ras, ras_class, none_left)

            addresses(i) = get_address(ras_class, string)

            symmetries(i) = get_abelian_sym(string)
        end do

        ! Check that the last string in the above loop really was the last string (as it should
        ! be) by calling the generating routine again.
        call generate_next_string(string, ras, ras_class, none_left)
        if (.not. none_left) call stop_all("find_symmetries", "Incorrect number of states found.")

        ! Now, sort the symmetries array from smallest to largest.
        call sort(symmetries, addresses)

        ! The addresses array was sorted with symmetries. If we have the final position of a state
        ! *after* the sort and we put that position into addresses, we get out the address of this
        ! state (as found by get_address). We want the inverse of this, a map that will take
        ! in an address from get_address and give out the position after the sort.
        do i = 1, ras_class%class_size
            ! g(f(x)) = x => g=f^(-1), so the following creates the inverse map.
            ras_class%address_map(addresses(i)) = i
        end do

        ! Count the number of states with each symmetry label.
        ras_class%num_sym = 0
        do i = 1, ras_class%class_size
            ras_class%num_sym(symmetries(i)) = ras_class%num_sym(symmetries(i)) + 1
        end do

        do i = 0, 7
            ras_class%cum_sym(i) = 0
            do j = 0, i - 1
                ras_class%cum_sym(i) = ras_class%cum_sym(i) + ras_class%num_sym(j)
            end do
        end do

        deallocate(symmetries)
        deallocate(addresses)

    end subroutine setup_ras_class