InitDavidsonCalc Subroutine

public subroutine InitDavidsonCalc(this, print_info, hamil_type)


Type IntentOptional Attributes Name
type(DavidsonCalcType), intent(inout) :: this
logical, intent(in) :: print_info
integer, intent(in) :: hamil_type


Source Code

Source Code

    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