subroutine find_and_output_lowdin_eigv(config_label, nvecs, overlap_matrix, hamil_matrix, npositive, &
all_evals, tOutput, spin_matrix, all_spin)
use CalcData, only: tWritePopsNorm, pops_norm
use util_mod, only: int_fmt, get_free_unit
integer, intent(in) :: config_label, nvecs
real(dp), intent(in) :: overlap_matrix(:, :), hamil_matrix(:, :)
integer, intent(out) :: npositive
real(dp), intent(out) :: all_evals(:, :)
logical, intent(in) :: tOutput
real(dp), optional, intent(in) :: spin_matrix(:, :)
real(dp), optional, intent(out) :: all_spin(:, :)
integer :: lwork, counter, i, nkeep, nkeep_len
integer :: info, ierr, temp_unit
real(dp), allocatable :: work(:)
real(dp), allocatable :: kp_final_hamil(:, :), kp_hamil_eigv(:), rotated_spin(:, :)
real(dp) :: kp_pert_energy_overlaps(nvecs)
character(7) :: string_fmt
character(25) :: ind1, filename
character(len=*), parameter :: stem = "lowdin"
if (tOutput) then
write(ind1, '(i15)') config_label
filename = trim(trim(stem)//'.'//trim(adjustl(ind1)))
temp_unit = get_free_unit()
open(temp_unit, file=trim(filename), status='replace')
if (tWritePopsNorm) then
write(temp_unit, '(4("-"),a41,25("-"))') "Norm of unperturbed initial wave function"
write(temp_unit, '(1x,es19.12,/)') sqrt(pops_norm)
end if
end if
! Create the workspace for the diagonaliser.
lwork = max(1, 3 * nvecs - 1)
allocate(work(lwork), stat=ierr)
kp_overlap_eigenvecs = overlap_matrix
! Now perform the diagonalisation.
call dsyev('V', 'U', nvecs, kp_overlap_eigenvecs, nvecs, kp_overlap_eigv, work, lwork, info)
npositive = 0
all_evals = 0.0_dp
if (present(spin_matrix)) all_spin = 0.0_dp
if (tOutput) write(temp_unit, '(4("-"),a26,40("-"))') "Overlap matrix eigenvalues"
do i = 1, nvecs
if (tOutput) write(temp_unit, '(1x,es19.12)') kp_overlap_eigv(i)
if (kp_overlap_eigv(i) > 0.0_dp) npositive = npositive + 1
end do
if (tOutput .and. tOverlapPert) then
write (temp_unit, '(/,"# The overlap of each final Hamiltonian eigenstate with each &
&requested perturbed ground state will be printed. The first &
&printed overlap is that with the first Krylov vector. The second &
&printed overlap is that with the vector specified with the &
&OVERLAP-PERTURB-ANNIHILATE and OVERLAP-PERTURB-CREATION &
&options.")')
end if
do nkeep = 1, npositive
allocate(kp_final_hamil(nkeep, nkeep))
allocate(kp_hamil_eigv(nkeep))
if (present(spin_matrix)) allocate(rotated_spin(nkeep, nkeep))
associate(transform_matrix => kp_transform_matrix(1:nvecs, 1:nkeep), &
inter_matrix => kp_inter_matrix(1:nvecs, 1:nkeep), &
eigenvecs_krylov => kp_eigenvecs_krylov(1:nvecs, 1:nkeep), &
init_overlaps => kp_init_overlaps(1:nkeep))
counter = 0
do i = nvecs - nkeep + 1, nvecs
counter = counter + 1
transform_matrix(:, counter) = kp_overlap_eigenvecs(:, i) / sqrt(kp_overlap_eigv(i))
end do
inter_matrix = matmul(hamil_matrix, transform_matrix)
kp_final_hamil = matmul(transpose(transform_matrix), inter_matrix)
call dsyev('V', 'U', nkeep, kp_final_hamil, nkeep, kp_hamil_eigv, work, lwork, info)
eigenvecs_krylov = matmul(transform_matrix, kp_final_hamil)
init_overlaps = matmul(overlap_matrix(1, :), eigenvecs_krylov) / scaling_factor
if (tOverlapPert) kp_pert_energy_overlaps(1:nkeep) = matmul(kp_all_pert_overlaps, eigenvecs_krylov)
! The spin matrix in the final eigenvalue basis.
if (present(spin_matrix)) then
! Use inter_matrix as temporary space.
inter_matrix = matmul(spin_matrix, eigenvecs_krylov)
rotated_spin = matmul(transpose(eigenvecs_krylov), inter_matrix)
end if
if (tOutput) then
nkeep_len = ceiling(log10(real(abs(nkeep) + 1, dp)))
write(string_fmt, '(i2,a5)') 15 - nkeep_len, '("-")'
write(temp_unit, '(/,4("-"),a37,'//int_fmt(nkeep, 1)//',1x,a12,'//string_fmt//')') &
"Eigenvalues and overlaps when keeping", nkeep, "eigenvectors"
do i = 1, nkeep
write(temp_unit, '(1x,es19.12,1x,es19.12)', advance='no') kp_hamil_eigv(i), init_overlaps(i)
if (tOverlapPert) write(temp_unit, '(1x,es19.12)', advance='no') kp_pert_energy_overlaps(i)
write(temp_unit, '()', advance='yes')
end do
end if
all_evals(1:nkeep, nkeep) = kp_hamil_eigv(1:nkeep)
if (present(spin_matrix)) then
do i = 1, nkeep
all_spin(i, nkeep) = rotated_spin(i, i)
end do
end if
end associate
deallocate(kp_final_hamil)
deallocate(kp_hamil_eigv)
if (allocated(rotated_spin)) deallocate(rotated_spin)
end do
deallocate(work)
if (tOutput) close(temp_unit)
end subroutine find_and_output_lowdin_eigv