subspace_expansion Subroutine

public subroutine subspace_expansion(this, basis_index)


Type IntentOptional Attributes Name
type(DavidsonCalcType), intent(inout) :: this
integer, intent(in) :: basis_index


Source Code

Source Code

    subroutine subspace_expansion(this, basis_index)
        type(DavidsonCalcType), intent(inout) :: this
        integer, intent(in) :: basis_index
        integer :: i
        real(dp) :: dot_prod, norm

        ! 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%super%space_size
            this%super%basis_vectors(i, basis_index) = this%residual(i) / (hamil_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
            dot_prod = dot_product(this%super%basis_vectors(:, basis_index), this%super%basis_vectors(:, i))
            this%super%basis_vectors(:, basis_index) = &
                this%super%basis_vectors(:, basis_index) - dot_prod * this%super%basis_vectors(:, i)
        end do

        ! 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%super%basis_vectors(:, basis_index), this%super%basis_vectors(:, basis_index))
        norm = sqrt(norm)
        this%super%basis_vectors(:, basis_index) = this%super%basis_vectors(:, basis_index) / norm

    end subroutine subspace_expansion