initialise_ras_space Subroutine

    subroutine initialise_ras_space(ras, classes)

        ! Given the five defining parameters for a RAS space, initialise the ras
        ! space by finding the ras classes corresponding to the ras space.

        type(ras_parameters), intent(inout) :: ras
        type(ras_class_data), intent(inout), allocatable, dimension(:) :: classes
        integer :: i, j, k, counter
        integer :: lower_ras3, upper_ras3

        tot_nelec = nel / 2
        tot_norbs = nbasis / 2

        ! Check that the RAS parameters are possible.
        if (ras%size_1 + ras%size_2 + ras%size_3 /= tot_norbs .or. &
            ras%min_1 > ras%size_1 * 2 .or. ras%max_3 > ras%size_3 * 2) &
            call stop_all("generate_ras", "RAS parameters are not possible.")
        if (mod(nel, 2) /= 0) call stop_all("generate_ras", "RAS-core only implmented for &
                                            & closed shell molecules.")

        ! First we need to find the different classes. A class is defined by the number
        ! of electrons in RAS1 and RAS3. Thus, we need to find all possible allowed
        ! combinations.

        ras%lower_ras1 = max(0, ras%min_1 - ras%size_1)
        ras%upper_ras1 = min(tot_nelec, ras%size_1)

        allocate(ras%class_label(ras%lower_ras1:ras%upper_ras1, 0:tot_nelec))
        ras%class_label = 0

        ! First count the number of RAS classes...
        counter = 0
        do i = ras%lower_ras1, ras%upper_ras1
            lower_ras3 = max(0, tot_nelec - i - ras%size_2)
            upper_ras3 = min(tot_nelec - i, ras%max_3)
            do j = lower_ras3, upper_ras3
                counter = counter + 1
                ras%class_label(i, j) = counter
            end do
        end do

        ras%num_classes = counter
        counter = 0

        ! ...then fill the classes in.
        allocate(classes(ras%num_classes))
        do i = ras%lower_ras1, ras%upper_ras1
            lower_ras3 = max(0, tot_nelec - i - 2 * ras%size_2)
            upper_ras3 = min(tot_nelec - i, ras%max_3)
            do j = lower_ras3, upper_ras3
                counter = counter + 1
                classes(counter)%nelec_1 = i
                classes(counter)%nelec_3 = j
                classes(counter)%nelec_2 = tot_nelec - i - j
                if (classes(counter)%nelec_2 < 0) then
                    call stop_all("initialise_ras_space", &
                                  "Current RAS combination is not possible.")
                end if
                allocate(classes(counter)%vertex_weights(0:tot_norbs, 0:tot_nelec))
                classes(counter)%vertex_weights = 0
            end do
        end do

        ! Form the vertex weights for each class.
        do i = 1, ras%num_classes
            classes(i)%vertex_weights(0, 0) = 1
            do j = 1, tot_norbs
                do k = 0, tot_nelec
                    ! If vertex not allowed then leave the corresponding vertex weight as 0.
                    if (vertex_not_allowed(classes(i)%nelec_1, &
                                           classes(i)%nelec_3, j, k, ras)) cycle
                    ! (Eq. 11.8.2)
                    if (k == 0) then
                        ! The first columns will always contain 1's.
                        classes(i)%vertex_weights(j, k) = 1
                    else
                        classes(i)%vertex_weights(j, k) = classes(i)%vertex_weights(j - 1, k) + &
                                                          classes(i)%vertex_weights(j - 1, k - 1)
                    end if
                end do
            end do
        end do

        ! Find the allowed combinations of classes in the full RAS space.
        do i = 1, ras%num_classes
            counter = 0
            allocate(classes(i)%allowed_combns(ras%num_classes))
            classes(i)%allowed_combns = 0
            do j = 1, ras%num_classes
                ! If the total number of electrons in the RAS spaces are correct with
                ! this combination.
                if (class_comb_allowed(ras, classes(i), classes(j))) then
                    counter = counter + 1
                    classes(i)%allowed_combns(counter) = j
                end if
            end do
            classes(i)%num_comb = counter
        end do

        allocate(ras%cum_classes(ras%num_classes))
        ras%cum_classes(1) = 0

        ras%num_strings = 0
        do i = 1, ras%num_classes
            call setup_ras_class(ras, classes(i))
            ras%num_strings = ras%num_strings + classes(i)%class_size
            if (i > 1) ras%cum_classes(i) = ras%cum_classes(i - 1) + classes(i - 1)%class_size
        end do

        HFSym_ras = int(HFSym%Sym%S)

    end subroutine initialise_ras_space