subspace_expansion_ss Subroutine

    subroutine subspace_expansion_ss(this, basis_index)

        type(davidson_ss), intent(inout) :: this
        integer, intent(in) :: basis_index
        integer :: i
        real(dp) :: dot_prod, dot_prod_tot, norm, norm_tot
        character(len=*), parameter :: t_r = "init_davidson_ss"

        ! Create the new basis state from the residual. This step performs
        ! t = (D - EI)^(-1) r,
        ! where D is the diagonal of the Hamiltonian matrix, E is the eigenvalue previously
        ! calculated, I is the identity matrix and r is the residual.

        do i = 1, this%space_size_this_proc
            this%basis_vectors(i, basis_index) = this%residual(i) / (cs_replicas(this%run)%core_ham_diag(i) - this%davidson_eigenvalue)
        end do

        ! This step then maskes the new basis vector orthogonal to all other basis vectors, by doing
        ! t <- t - (t,v)v
        ! for each basis vector v, where (t,v) denotes the dot product.
        do i = 1, basis_index - 1
            if (this%space_size_this_proc > 0) then
                dot_prod = dot_product(this%basis_vectors(:, basis_index), this%basis_vectors(:, i))
            else
                dot_prod = 0.0_dp
            end if

            call MPISumAll(dot_prod, dot_prod_tot)

            if (this%space_size_this_proc > 0) then
                this%basis_vectors(:, basis_index) = this%basis_vectors(:, basis_index) - dot_prod_tot * this%basis_vectors(:, i)
            end if
        end do

        if (this%space_size_this_proc > 0) then
            ! Finally we calculate the norm of the new basis vector and then normalise it to have a norm of 1.
            ! The new basis vector is stored in the next available column in the basis_vectors array.
            norm = dot_product(this%basis_vectors(:, basis_index), this%basis_vectors(:, basis_index))
        else
            norm = 0.0_dp
        end if

        call MPISumAll(norm, norm_tot)
        norm = sqrt(norm_tot)

        if (this%space_size_this_proc > 0) then
            this%basis_vectors(:, basis_index) = this%basis_vectors(:, basis_index) / norm
        end if

    end subroutine subspace_expansion_ss