subroutine init_davidson_ss(this, print_info, run)
! 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 util_mod, only: int_fmt
type(davidson_ss), intent(inout) :: this
logical, intent(in) :: print_info
integer, intent(in) :: run
integer :: i, hfindex, hf_proc, mem_reqd, mem_reqd_full, ierr
real(dp) :: hf_elem, hf_elem_this_proc, hf_elem_all_procs(0:nProcessors - 1)
logical :: skip_calc, skip_calc_all(0:nProcessors - 1)
integer(MPIArg) :: mpi_temp
character(len=*), parameter :: t_r = "init_davidson_ss"
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
this%run = run
associate ( &
davidson_eigenvalue => this%davidson_eigenvalue, &
space_size => this%space_size, &
space_size_this_proc => this%space_size_this_proc, &
rep => cs_replicas(this%run) &
)
space_size_this_proc = size(rep%core_ham_diag)
allocate(this%displs(0:nProcessors - 1))
allocate(this%sizes(0:nProcessors - 1))
mpi_temp = int(space_size_this_proc, MPIArg)
call MPIAllGather(mpi_temp, this%sizes, ierr)
! The total space size across all processors.
space_size = int(sum(this%sizes))
this%displs(0) = 0
do i = 1, nProcessors - 1
!this%displs(i) = sum(this%displs(:i-1))
this%displs(i) = this%displs(i - 1) + this%sizes(i - 1)
end do
! 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)
end if
safe_calloc_e(this%davidson_eigenvector, (space_size_this_proc), 0.0_dp, 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_temp(1)
! call mpibcast(davidson_eigenvalue)
! skip_calc = .true.
! return
!end if
if (print_info) then
write(stdout,*) 'Space sizes and max Davidson iterations: ', space_size_this_proc, max_num_davidson_iters
call neci_flush(stdout)
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_this_proc * 8) / 1000000
! the memory required to allocate residual.
mem_reqd_full = space_size * 8 / 1000000
if (print_info) then
write(stdout, '(1x,"allocating array to hold subspace vectors (",'//int_fmt(mem_reqd, 0)//',1x,"mb).")') mem_reqd
call neci_flush(stdout)
end if
if (print_info) then
write (stdout, '(1x,"allocating array to hold multiplied krylov vectors (",' &
//int_fmt(mem_reqd, 0)//',1x,"mb).")') mem_reqd
call neci_flush(stdout)
end if
if (print_info) then
write (stdout, '(1x,"allocating temporary vector (",' &
//int_fmt(mem_reqd_full, 0)//',1x,"mb).",/)') mem_reqd_full
call neci_flush(stdout)
end if
safe_calloc(this%projected_hamil, (max_num_davidson_iters, max_num_davidson_iters), 0.0_dp)
safe_calloc(this%projected_hamil_work, (max_num_davidson_iters, max_num_davidson_iters), 0.0_dp)
hf_elem_this_proc = maxval(-rep%core_ham_diag)
call MPIAllGather(hf_elem_this_proc, hf_elem_all_procs, ierr)
! Find the processor on which the HF determinant lives:
hf_proc = maxloc((hf_elem_all_procs), 1) - 1
! The Hartree--Fock element itself
hf_elem = -maxval(hf_elem_all_procs)
! allocate the necessary arrays:
safe_calloc(this%basis_vectors, (space_size_this_proc, max_num_davidson_iters), 0.0_dp)
safe_calloc(this%multiplied_basis_vectors, (space_size_this_proc, max_num_davidson_iters), 0.0_dp)
safe_calloc(this%residual, (space_size_this_proc), 0.0_dp)
safe_calloc(this%eigenvector_proj, (max_num_davidson_iters), 0.0_dp)
safe_calloc(this%full_vector, (space_size), 0.0_dp)
! If the HF determinant is on this process:
if (hf_proc == iProcIndex) then
hfindex = maxloc((-rep%core_ham_diag), 1)
! for the initial basis vector, choose the hartree-fock state:
!this%super%basis_vectors(hfindex, 1) = 1.0_dp
this%basis_vectors(hfindex, 1) = 1.0_dp
! choose the hartree-fock state as the initial guess at the ground state, too.
this%davidson_eigenvector(hfindex) = 1.0_dp
end if
! Set the initial eigenvector in the Davidson basis - there's only one
! Davidson vector to start with, so this is trivial.
this%eigenvector_proj(1) = 1.0_dp
! fill in the projected hamiltonian so far.
this%projected_hamil(1, 1) = hf_elem
! 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 = hf_elem - 0.001_dp
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.
call multiply_hamil_and_vector_ss(this%davidson_eigenvector, this%multiplied_basis_vectors(:, 1), &
this%full_vector, this%sizes, this%displs, &
this%run)
if (space_size_this_proc > 0) then
if (all(abs(this%multiplied_basis_vectors(:, 1) - hf_elem * this%davidson_eigenvector) < 1.0e-12_dp)) then
skip_calc = .true.
end if
else
skip_calc = .true.
end if
call MPIAllGather(skip_calc, skip_calc_all, ierr)
if (all(skip_calc_all)) return
! calculate the intial residual vector.
call calculate_residual_ss(this, 1)
call calculate_residual_norm_ss(this)
if (print_info) write(stdout, '(1x,"done.",/)'); call neci_flush(stdout)
end associate
end subroutine init_davidson_ss