subroutine perform_lanczos(this, det_list, n_states, hamil_type, print_info_in)
use CalcData, only: t_lanczos_orthogonalise, t_lanczos_store_vecs, &
lanczos_max_restarts, lanczos_max_vecs, &
lanczos_energy_precision, lanczos_ritz_overlap_precision
use hamiltonian_linalg, only: orthogonalise_against_previous_basis_vectors
type(LanczosCalcType), intent(inout) :: this
integer, intent(in) :: det_list(:, :), n_states, hamil_type
logical, intent(in) :: print_info_in
logical :: print_info, t_deltas_pass
integer :: k, i, ierr, delta_check_outcome, restart_counter
real(dp) :: start_time, end_time
real(dp) :: exp_val, overlap
character(40) :: main_output_fmt, final_output_fmt
character(*), parameter :: t_r = "perform_lanczos"
if (.not. t_lanczos_store_vecs .and. t_lanczos_orthogonalise) then
call stop_all(t_r, "storage of Lanczos method basis vectors needed for orthogonalisation")
end if
! format specifier formatting!
! ensure that enough places are displayed to show convergence to the desired accuracy
write(main_output_fmt, '("(8X,i4,3X,i2,2x,f",i2,".",i2,",2x,f9.3)")') lanczos_energy_precision + 7, lanczos_energy_precision
write(final_output_fmt, '("(1x,",A7,",i2,4X,f",i2,".",i2,")")') '"State"', lanczos_energy_precision + 7, lanczos_energy_precision
! Only let the root processor print information.
print_info = print_info_in .and. (iProcIndex == root)
if (.not. allocated(this%lanczos_vector)) then
call InitLanczosCalc(this, det_list, print_info, hamil_type, n_states, &
lanczos_max_vecs, t_lanczos_store_vecs, t_lanczos_orthogonalise, &
lanczos_max_restarts, lanczos_energy_precision, lanczos_ritz_overlap_precision)
end if
if (print_info) then
write(stdout, '(1X,"Perfoming a Lanczos Diagonalisation of the trial space")')
if (this%super%t_orthogonalise) then
write(stdout, '(/,1X,"Orthogonalising Lanczos vectors")')
else
write(stdout, '(/,1X,"Not orthogonalising Lanczos vectors")')
end if
end if
associate ( &
basis_vectors => this%super%basis_vectors, &
current_v => this%current_v, &
old_v => this%old_v &
)
! start Lanczos restart loop
restart_counter = 0
do
if (iprocindex == root) then
call addVec(this, 1, this%lanczos_vector)
end if
call MPIBCast(this%lanczos_vector)
call MPIBCast(this%super%basis_vectors(:, 1))
! set lanczos vector 1 to H*first_v
call project_hamiltonian_lanczos(this, 1)
if (print_info) then
write(stdout, '(/,1X,"Iteration",4x,"State",12X,"Energy",7X,"Time")'); call neci_flush(stdout)
end if
! start Lanczos main loop
do k = 1, this%super%max_subspace_size - 1
if (this%super%skip_calc) exit
start_time = MPI_WTIME()
if (iprocindex == root) then
if (this%super%t_orthogonalise .and. t_lanczos_store_vecs) then
! Although Lanczos vectors computed in exact arithmetic are orthogonal,
! those generated on a machine will stray increasingly from orthogonal with
! each iteration: so do GS procedure here
call orthogonalise_against_previous_basis_vectors(this%super, k + 1)
end if
if (t_lanczos_store_vecs) then
overlap = real(inner_product(basis_vectors(:, k), this%lanczos_vector), dp)
call setAlpha(this, k, overlap)
call addVec(this, k + 1, this%lanczos_vector - getAlpha(this, k) * basis_vectors(:, k))
call setBeta(this, k + 1, euclidean_norm(basis_vectors(:, k + 1)))
basis_vectors(:, k + 1) = basis_vectors(:, k + 1) / getBeta(this, k + 1)
else
overlap = real(dot_product(this%current_v, this%lanczos_vector), dp)
call setAlpha(this, k, overlap)
call addVec(this, k + 1, this%lanczos_vector - getAlpha(this, k) * current_v)
call setBeta(this, k + 1, sqrt(real(dot_product(current_v, current_v), dp)))
current_v = current_v / getBeta(this, k + 1)
end if
end if
if (t_lanczos_store_vecs) then
call MPIBCast(this%super%basis_vectors(:, k + 1))
else
call MPIBCast(this%current_v)
end if
call project_hamiltonian_lanczos(this, k + 1)
end_time = MPI_WTIME()
if (iprocindex == root) then
if (t_lanczos_store_vecs) then
this%lanczos_vector = this%lanczos_vector - getBeta(this, k + 1) * basis_vectors(:, k)
else
this%lanczos_vector = this%lanczos_vector - getBeta(this, k + 1) * old_v
end if
if (t_non_hermitian_2_body) then
call diagonalise_tridiagonal_non_hermitian(this, k, .false.)
else
call diagonalise_tridiagonal(this, k, .false.)
end if
! if all states have converged we can end the main loop prematurely
t_deltas_pass = .true.
do i = 1, this%n_states
t_deltas_pass = t_deltas_pass .and. check_delta(this, k, i)
if (.not. this%t_states_converged(i)) then
if (print_info) then
write(stdout, trim(main_output_fmt)) k, i, this%ritz_values(i), end_time - start_time
call neci_flush(stdout)
end if
end if
end do
write(stdout, *)
end if
call MPIBCast(t_deltas_pass)
if (t_deltas_pass) then
exit
end if
call MPIBCast(this%lanczos_vector)
! end Lanczos main loop
end do
if (iprocindex == root) then
if (t_non_hermitian_2_body) then
call diagonalise_tridiagonal_non_hermitian(this, k, .true.)
else
call diagonalise_tridiagonal(this, k, .true.)
end if
call compute_ritz_vectors(this, k)
end if
! we now evaluate the quality of the convergence achieved. If a state is sufficiently
! well converged, save it and its energy otherwise, add the current best ritz vector
! for the state to the restarting Lanczos vector such that we restart with an equal
! superposition of the best guesses for the unconverged states.
if (iProcIndex == root) then
this%lanczos_vector(:) = h_cast(0.0_dp)
do i = 1, this%n_states
if (this%t_states_converged(i)) then
cycle
end if
if (.not. check_delta(this, k, i)) then
this%lanczos_vector = this%lanczos_vector + this%ritz_vectors(:, i)
else
this%eigenvectors(:, i) = this%ritz_vectors(:, i)
this%eigenvalues(i) = this%ritz_values(i)
this%t_states_converged(i) = .true.
end if
end do
end if
call MPIBCast(this%t_states_converged)
if (.not. any(.not. this%t_states_converged)) then
exit
end if
! only normalize the lanczos vector if we are not
! already converged, as it will be 0 else
this%lanczos_vector = this%lanczos_vector / euclidean_norm(this%lanczos_vector)
restart_counter = restart_counter + 1
if (restart_counter < this%max_restarts) then
if (print_info) then
write(stdout, '(/,1x,"Maximum iteration number reached, restarting Lanczos algorithm (",i3,"/",i3")")') &
restart_counter, this%max_restarts
end if
else
if (print_info) then
write(stdout, '(/,1x,"Maximum restart number reached. Some states may not be converged.")')
end if
exit
end if
! end Lanczos restart loop
end do
! serialize the printout
if (iprocindex == root) then
do i = 1, this%n_states
this%eigenvectors(:, i) = this%ritz_vectors(:, i)
this%eigenvalues(i) = this%ritz_values(i)
end do
if (check_deltas(this, k)) then
if (print_info) then
write(stdout, '(i2" eigenvalues(s) were successfully converged to within ",5ES16.7)') &
this%n_states, this%convergence_error
call neci_flush(stdout)
end if
end if
call perform_orthogonality_test(this)
if (print_info) then
write(stdout, '(/,1x,"Final calculated energies:")')
do i = 1, this%n_states
write(stdout, final_output_fmt) i, this%eigenvalues(i)
call neci_flush(stdout)
end do
end if
end if
! how good are the ritz vectors?
write(stdout, '(/,1x,"Ritz vector expectation energies:")')
do i = 1, this%n_states
exp_val = get_rayleigh_quotient(this, i)
if (print_info) then
write(stdout, final_output_fmt) i, exp_val
call neci_flush(stdout)
end if
end do
write(stdout, '(/,1x,"Ritz vector residual norms:")')
do i = 1, this%n_states
exp_val = compute_residual_norm(this, i)
if (print_info) then
write(stdout, final_output_fmt) i, exp_val
call neci_flush(stdout)
end if
end do
write(stdout, '(/,1x,"End of Lanczos procedure.",/)')
! wait for final diagonalisation before freeing any memory
call MPIBarrier(ierr)
call FreeLanczosCalc(this)
end associate
end subroutine perform_lanczos