davidson_neci.F90 Source File


Source Code

Source Code

#include "macros.h"
module davidson_neci

    ! This module performs the Davidson method to find the ground state of a diagonally-
    ! dominant matrix. For details of the theory behind the method, see i.e:
    ! http://web.mit.edu/bolin/www/Project-Report-18.335J.pdf

    use constants
    use FciMCData, only: hamiltonian, DavidsonTag, user_input_max_davidson_iters, &
    use MemoryManager, only: TagIntType, LogMemAlloc, LogMemDealloc
    use Parallel_neci, only: iProcIndex, nProcessors, MPIArg, MPIBarrier
    use Parallel_neci, only: MPIBCast, MPIGatherV, MPIAllGather
    use MPI_wrapper, only: root
    use ras_data
    use sparse_arrays, only: sparse_ham, hamil_diag, HDiagTag
    use util_mod, only: neci_flush, stop_all
    use hamiltonian_linalg, only: &
        full_hamil_type, &
        sparse_hamil_type, &
        parallel_sparse_hamil_type, &
        direct_ci_type, &
        HamiltonianCalcType, &
        initHamiltonianCalc, &
        multiply_hamil_and_vector, &
        direct_ci_inp, &
        direct_ci_out, &

    implicit none

    integer :: max_num_davidson_iters = 50
    real(dp) :: residual_norm_target = 1e-7_dp

    ! To cut down on the amount of global data, introduce a derived type to hold a Davidson session
    type DavidsonCalcType
        ! "super type"
        type(HamiltonianCalcType) :: super
        ! This array stores the basis vectors multiplied by H in its columns, i.e.
        ! multiplied_basis_vectors(:,1) = H*basis_vector(:,1).
        real(dp), allocatable, dimension(:, :) :: multiplied_basis_vectors
        ! By diagonalising the projected Hamiltonian we get an estimate at the ground state in
        ! the basis of those basis vectors stored in the basis_vectors array. davidson_eigenvector
        ! stores this same state, but in the *original* basis set. It therefore has a dimension
        ! the same size as the vector space.
        real(dp), allocatable, dimension(:) :: davidson_eigenvector
        ! This array holds the components of davidson_eigenvector in the basis of Krylov vectors.
        real(dp), allocatable, dimension(:) :: eigenvector_proj
        ! The residual is defined as r = H*v - E*v, where H is the Hamiltonian matrix, v is the
        ! ground state estimate (stored in davidson_eigenvector) and E is the corresponding
        ! energy eigenvalue. If v is an exact eigenstate then all the components of the residual
        ! are zero.
        real(dp), allocatable, dimension(:) :: residual
        ! As noted above, if davidson_eigenvector holds an exact eigenstate then the residual
        ! will have all zero components and this norm (the standard Euclidean norm) will be zero.
        ! Hence it is a measure of how converged the solution is.
        real(dp) :: residual_norm
        real(dp) :: davidson_eigenvalue
        ! temp vectors for real matrix-vector calculations even when compiling in complex mode
        real(dp), allocatable, dimension(:) :: temp_in, temp_out
    end type


    subroutine perform_davidson(this, hamil_type_in, print_info_in)
        use mpi
        integer, intent(in) :: hamil_type_in
        logical, intent(in) :: print_info_in
        logical :: print_info
        integer :: i
        real(dp) :: start_time, end_time
        type(DavidsonCalcType), intent(inout) :: this
        character(*), parameter :: this_routine = "perform_davidson"

        ! Only let the root processor print information.
        print_info = print_info_in .and. (iProcIndex == root)

        call InitDavidsonCalc(this, print_info, hamil_type_in)

        if (print_info) write(stdout, '(1X,"Iteration",4X,"Residual norm",12X,"Energy",7X,"Time")'); call neci_flush(stdout)

        do i = 2, max_num_davidson_iters

            if (this%super%skip_calc) exit

            start_time = MPI_WTIME()

            if (iProcIndex == root) call subspace_expansion(this, i)

            call project_hamiltonian(this, i)

            if (iProcIndex == root) call subspace_extraction(this, i)

            call calculate_residual(this, i)

            call calculate_residual_norm(this)

            end_time = MPI_WTIME()

            if (print_info) write(stdout, '(8X,i2,3X,f14.9,2x,f16.10,2x,f9.3)') i - 1, this%residual_norm, &
                this%davidson_eigenvalue, end_time - start_time; call neci_flush(stdout)

            if (this%residual_norm < residual_norm_target) exit

            if (i == max_num_davidson_iters) then
                call stop_all(this_routine, "Davidson iteration reached the maximum number of iterations. &
                                            &The deterministic energy may not be converged. &
                                            &You can increase 'davidson-max-iters' or 'davidson-target-tolerance'.")
            end if

        end do

        if (print_info) write(stdout, '(/,1x,"Final calculated correlation energy:",1X,f16.10)') this%davidson_eigenvalue

        call FreeDavidsonCalc(this)

    end subroutine perform_davidson

    subroutine InitDavidsonCalc(this, print_info, hamil_type)

        ! This subroutine initialises the Davdison method by allocating the necessary arrays,
        ! defining the initial basis vector and projected Hamiltonian, and setting an initial
        ! guess at the ground state eigenvalue. It also calculates the corresponding residual
        ! which is needed to expand the space.

        use direct_ci, only: create_ham_diag_direct_ci
        use FciMCData, only: davidson_ras, davidson_classes, davidson_strings
        use ras, only: find_ras_size
        use util_mod, only: int_fmt

        type(DavidsonCalcType), intent(inout) :: this

        logical, intent(in) :: print_info
        integer, intent(in) :: hamil_type

        logical :: skip_calc
        integer :: i, mem_reqd, residual_mem_reqd, ierr
        integer(mpiarg) :: mpi_temp
        real(dp), allocatable :: hamil_diag_temp(:)
        character(len=*), parameter :: t_r = "init_davidson"

        if( allocated(user_input_max_davidson_iters) ) &
            max_num_davidson_iters = user_input_max_davidson_iters

        if( allocated(user_input_davidson_tolerance) ) &
            residual_norm_target = user_input_davidson_tolerance

        call InitHamiltonianCalc(this%super, print_info, hamil_type, max_num_davidson_iters, .true., .false.)

        associate ( &
            davidson_eigenvalue => this%davidson_eigenvalue, &
            space_size => this%super%space_size, &
            hfindex => this%super%hfindex &

            ! if a davidson calculation has already been performed, this array might still be
            ! allocated, so check!
            if (allocated(this%davidson_eigenvector)) then
                deallocate(this%davidson_eigenvector, stat=ierr)
                call logmemdealloc(t_r, davidsontag, ierr)
            end if
            safe_calloc_e(this%davidson_eigenvector, (space_size), 0.0_dp, ierr)
            call logmemalloc("davidson_eigenvector", space_size, 8, t_r, davidsontag, ierr)

            ! if there is only one state in the space being diagonalised:
            if (space_size == 1) then
                this%davidson_eigenvector(1) = 1.0_dp
                if (iprocindex == root) davidson_eigenvalue = hamil_diag(1)
                call mpibcast(davidson_eigenvalue)
                skip_calc = .true.
                call mpibcast(skip_calc)
                this%super%skip_calc = skip_calc
            end if

            if (iprocindex == root) then
                if (tCalcHFIndex) then
                    hfindex = maxloc((-hamil_diag), 1)
                end if
                ! the memory required to allocate each of basis_vectors and
                ! multipied_basis_vectors, in mb.
                mem_reqd = (max_num_davidson_iters * space_size * 8) / 1000000
                ! the memory required to allocate residual.
                residual_mem_reqd = space_size * 8 / 1000000

                ! allocate the necessary arrays:
                if (print_info) write (stdout, '(1x,"allocating array to hold multiplied krylov vectors (",' &
                                       //int_fmt(mem_reqd, 0)//',1x,"mb).")') mem_reqd; call neci_flush(stdout)
                safe_calloc(this%multiplied_basis_vectors, (space_size, max_num_davidson_iters), 0.0_dp)
                safe_calloc(this%eigenvector_proj, (max_num_davidson_iters), 0.0_dp)
                if (print_info) write (stdout, '(1x,"allocating array to hold the residual vector (",' &
                                       //int_fmt(residual_mem_reqd, 0)//',1x,"mb).",/)') residual_mem_reqd; call neci_flush(stdout)
                safe_calloc(this%residual, (space_size), 0.0_dp)

                ! for the initial basis vector, choose the hartree-fock state:
                this%super%basis_vectors(hfindex, 1) = 1.0_dp
                ! choose the hartree-fock state as the initial guess at the ground state, too.
                this%eigenvector_proj(1) = 1.0_dp
                this%davidson_eigenvector(hfindex) = 1.0_dp

                ! fill in the projected hamiltonian so far.
                this%super%projected_hamil(1, 1) = hamil_diag(hfindex)
                ! take the initial eigenvalue to be the hartree-fock energy minus some small
                ! amount. this value cannot be exactly the hartree-fock energy, as this will
                ! result in dividing by zero in the subspace expansion step.
                davidson_eigenvalue = hamil_diag(hfindex) - 0.001_dp
                safe_malloc(this%temp_in, (space_size))
                safe_malloc(this%temp_out, (space_size))
            end if

            if (print_info) write(stdout, '(1x,"calculating the initial residual vector...")', advance='no'); call neci_flush(stdout)

            ! check that multiplying the initial vector by the hamiltonian doesn't give back
            ! the same vector. if it does then the initial vector (the hf determinant) is
            ! the ground state, so just keep that and exit the calculation.
            ! also, the result of the multiplied basis vector is used to calculate the
            ! initial residual vector, if the above condition is not true.
            skip_calc = .false.
            if (iprocindex == root) then
                call multiply_hamil_and_vector(this%super, this%davidson_eigenvector, this%multiplied_basis_vectors(:, 1))
                call multiply_hamil_and_vector(this%super, this%davidson_eigenvector, this%temp_out)
            end if
            if (iprocindex == root) then
                if (all(abs(this%multiplied_basis_vectors(:, 1) - hamil_diag(hfindex) * this%davidson_eigenvector) < 1.0e-12_dp)) then
                    skip_calc = .true.
                    davidson_eigenvalue = hamil_diag(hfindex)
                end if
            end if
            if (hamil_type == parallel_sparse_hamil_type) call mpibcast(skip_calc)
            this%super%skip_calc = skip_calc
            if (skip_calc) return
            ! calculate the intial residual vector.
            call calculate_residual(this, 1)
            call calculate_residual_norm(this)

            if (print_info) write(stdout, '(1x,"done.",/)'); call neci_flush(stdout)

        end associate

    end subroutine InitDavidsonCalc

    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

    subroutine subspace_extraction(this, basis_index)

        type(DavidsonCalcType), intent(inout) :: this
        integer, intent(in) :: basis_index
        integer :: lwork, info
        real(dp), allocatable, dimension(:) :: work
        real(dp) :: eigenvalue_list(basis_index)
        ! Scrap space for the diagonaliser.
        lwork = max(1, 3 * basis_index - 1)

        ! This routine diagonalises a symmetric matrix, A.
        ! V tells the routine to calculate eigenvalues *and* eigenvectors.
        ! U tells the routine to get the upper half of A (it is symmetric).
        ! basis_index is the number of rows and columns in A.
        ! A = projected_hamil_work. This matrix stores the eigenvectors in its columns on output.
        ! basis_index is the leading dimension of A.
        ! eigenvalue_list stores the eigenvalues on output.
        ! work is scrap space.
        ! lwork is the length of the work array.
        ! info = 0 on output is diagonalisation is successful.
        call dsyev( &
            'V', &
            'U', &
            basis_index, &
            this%super%projected_hamil_work(1:basis_index, 1:basis_index), &
            basis_index, &
            eigenvalue_list, &
            work, &
            lwork, &
            info &

        ! The first column stores the ground state.
        this%davidson_eigenvalue = eigenvalue_list(1)
        this%eigenvector_proj(1:basis_index) = this%super%projected_hamil_work(1:basis_index, 1)

        ! eigenvector_proj stores the eigenstate in the basis of vectors stored in the array
        ! basis_vectors. We now want it in terms of the original basis. To get this, multiply
        ! eigenvector_proj by basis_vectors(:, 1:basis_index):

        ! This function performs y := alpha*a*x + beta*y
        ! N specifies not to use the transpose of a.
        ! space_size is the number of rows in a.
        ! basis_index is the number of columns of a.
        ! alpha = 1.0_dp.
        ! a = basis_vectors(:,1:basis_index).
        ! space_size is the first dimension of a.
        ! input x = eigenvector_proj(1:basis_index).
        ! 1 is the increment of the elements of x.
        ! beta = 0.0_dp.
        ! output y = davidson_eigenvector.
        ! 1 is the incremenet of the elements of y.
        call dgemv('N', &
                   this%super%space_size, &
                   basis_index, &
                   1.0_dp, &
                   this%super%basis_vectors(:, 1:basis_index), &
                   this%super%space_size, &
                   this%eigenvector_proj(1:basis_index), &
                   1, &
                   0.0_dp, &
                   this%davidson_eigenvector, &

    end subroutine subspace_extraction

    subroutine project_hamiltonian(this, basis_index)

        type(DavidsonCalcType), intent(inout) :: this
        integer, intent(in) :: basis_index
        integer :: i

        if (iProcIndex == root) then
            ! Multiply the new basis_vector by the hamiltonian and store the result in
            ! multiplied_basis_vectors.
            call multiply_hamil_and_vector(this%super, real(this%super%basis_vectors(:, basis_index), dp), &
                                           this%multiplied_basis_vectors(:, basis_index))

            ! Now multiply U^T by (H U) to find projected_hamil. The projected Hamiltonian will
            ! only differ in the new final column and row. Also, projected_hamil is symmetric.
            ! Hence, we only need to calculate the final column, and use this to update the final
            ! row also.
            do i = 1, basis_index
                this%super%projected_hamil(i, basis_index) = &
                    dot_product(this%super%basis_vectors(:, i), this%multiplied_basis_vectors(:, basis_index))
                this%super%projected_hamil(basis_index, i) = this%super%projected_hamil(i, basis_index)
            end do

            ! We will use the scrap Hamiltonian to pass into the diagonaliser later, since it
            ! overwrites this input matrix with the eigenvectors. Hence, make sure the scrap space
            ! stores the updated projected Hamiltonian.
            this%super%projected_hamil_work = this%super%projected_hamil
            call multiply_hamil_and_vector(this%super, this%temp_in, this%temp_out)
        end if

    end subroutine project_hamiltonian

    subroutine calculate_residual(this, basis_index)

        type(DavidsonCalcType), intent(inout) :: this
        integer, intent(in) :: basis_index

        ! This routine calculates the residual, r, corresponding to the new estimate of the
        ! ground state, stored in davidson_eigenvector. This is defined as
        ! r = Hv - Ev,
        ! where H is the Hamiltonian, v is the ground state vector estimate and E is the
        ! ground state energy estimate.

        if (iProcIndex == root) then
            ! Calculate r = Hv - Ev:
            ! Note that, here, eigenvector_proj holds the components of v in the Krylov basis,
            ! and multiplied_basis_vectors holds the Krylov vectors multiplied by H, hence
            ! the matmul below does indeed retturn Hv.
            this%residual = matmul(this%multiplied_basis_vectors(:, 1:basis_index), this%eigenvector_proj(1:basis_index))
            this%residual = this%residual - this%davidson_eigenvalue * this%davidson_eigenvector
        end if

    end subroutine calculate_residual

    subroutine calculate_residual_norm(this)

        type(DavidsonCalcType), intent(inout) :: this
        ! This subroutine calculates the Euclidean norm of the reisudal vector, r:
        ! residual_norm^2 = \sum_i r_i^2
        if (iProcIndex == root) then
            this%residual_norm = sqrt(dot_product(this%residual, this%residual))
        end if

        if (this%super%hamil_type == parallel_sparse_hamil_type) call MPIBCast(this%residual_norm)

    end subroutine calculate_residual_norm

    subroutine FreeDavidsonCalc(this)
        use hamiltonian_linalg, only: DestroyHamiltonianCalc
        type(DavidsonCalcType), intent(inout) :: this
        ! destroy the super type instance
        call DestroyHamiltonianCalc(this%super)
        ! we are now done with these arrays
        ! but keep the davidson eigenvector
    end subroutine FreeDavidsonCalc

    subroutine DestroyDavidsonCalc(this)
        type(DavidsonCalcType), intent(inout) :: this
        ! deallocate the davidson vector as well
        call FreeDavidsonCalc(this)
    end subroutine DestroyDavidsonCalc

    function davidson_direct_ci_init() result(this)
        use bit_rep_data, only: NIfD
        use direct_ci, only: create_direct_ci_arrays
        use FciMCData, only: davidson_ras, davidson_classes, davidson_strings, davidson_iluts, davidson_excits
        use ras, only: initialise_ras_space, find_ras_size

        type(DavidsonCalcType) :: this
        integer :: class_i, class_j, j, sym_i, sym_j

        write(stdout, '(/,1X,"Beginning Direct CI Davidson calculation.",/)'); call neci_flush(stdout)

        call initialise_ras_space(davidson_ras, davidson_classes)
        ! The total hilbert space dimension of calculation to be performed.
        call find_ras_size(davidson_ras, davidson_classes, this%super%space_size)

        allocate(davidson_strings(-1:tot_nelec, davidson_ras%num_strings))
        allocate(davidson_iluts(0:NIfD, davidson_ras%num_strings))

        ! Create the arrays used by the direct CI multiplication.
        call create_direct_ci_arrays(davidson_ras, davidson_classes, davidson_strings, &
                                     davidson_iluts, davidson_excits)

        ! Allocate input and output direct CI vectors.
        allocate(direct_ci_inp(size(davidson_classes), size(davidson_classes), 0:7))
        allocate(direct_ci_out(size(davidson_classes), size(davidson_classes), 0:7))
        do class_i = 1, size(davidson_classes)
            do j = 1, davidson_classes(class_i)%num_comb
                class_j = davidson_classes(class_i)%allowed_combns(j)
                do sym_i = 0, 7
                    sym_j = ieor(int(HFSym_ras), sym_i)
                    if (davidson_classes(class_i)%num_sym(sym_i) == 0) cycle
                    if (davidson_classes(class_j)%num_sym(sym_j) == 0) cycle
                    allocate(direct_ci_inp(class_i, class_j, sym_i)% &
                              elements(1:davidson_classes(class_i)%num_sym(sym_i), 1:davidson_classes(class_j)%num_sym(sym_j)))
                    allocate(direct_ci_out(class_i, class_j, sym_i)% &
                              elements(1:davidson_classes(class_i)%num_sym(sym_i), 1:davidson_classes(class_j)%num_sym(sym_j)))
                end do
            end do
        end do
    end function davidson_direct_ci_init

    subroutine davidson_direct_ci_end(this)
        type(DavidsonCalcType), intent(inout) :: this
        integer :: ierr
        if (allocated(hamil_diag)) then
            deallocate(hamil_diag, stat=ierr)
            call LogMemDealloc("davidson_direct_ci_end", HDiagTag, ierr)
        end if
        if (allocated(this%davidson_eigenvector)) then
            deallocate(this%davidson_eigenvector, stat=ierr)
            call LogMemDealloc("davidson_direct_ci_end", DavidsonTag, ierr)
        end if

        write(stdout, '(/,1X,"Direct CI Davidson calculation complete.",/)'); call neci_flush(stdout)

        write(stdout, "(1X,a10,f16.10)") "GROUND E =", this%davidson_eigenvalue; call neci_flush(stdout)
    end subroutine davidson_direct_ci_end

end module davidson_neci