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
return
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
else
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))
else
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