subroutine find_and_output_gs_eigv(config_label, nvecs)
! Use the Gram-Schmidt approach rather than the Lowdin approach above
! (see Phys. Rev. B. 85, 205119).
use util_mod, only: int_fmt, get_free_unit
integer, intent(in) :: config_label, nvecs
integer :: lwork, counter, nkeep, nkeep_len, temp_unit
integer :: npositive, info, ierr
integer :: i, j, n, m
real(dp), allocatable :: work(:)
real(dp), allocatable :: kp_final_hamil(:, :), kp_hamil_eigv(:)
character(7) :: string_fmt
character(25) :: ind1, filename
character(len=*), parameter :: stem = "gram_schmidt"
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')
! Create the workspace for the diagonaliser.
lwork = max(1, 3 * nvecs - 1)
allocate(work(lwork), stat=ierr)
! Use the following allocated arrays as work space for the following routine. Not ideal, I know...
allocate(kp_final_hamil(nvecs, nvecs))
associate(S_tilde => kp_inter_matrix, N => kp_init_overlaps)
call construct_gs_transform_matrix(kp_overlap_mean, kp_transform_matrix, S_tilde, kp_final_hamil, N, nvecs)
npositive = 0
do i = 1, nvecs
if (N(i) > 0.0_dp) then
npositive = npositive + 1
else
exit
end if
end do
write(temp_unit, '(4("-"),a21,45("-"))') "Normalisation factors"
do i = 1, npositive
write(temp_unit, '(1x,es19.12)') N(i)
end do
end associate
deallocate(kp_final_hamil)
do nkeep = 1, npositive
allocate(kp_final_hamil(nkeep, nkeep))
allocate(kp_hamil_eigv(nkeep))
associate(S => kp_transform_matrix(1:nvecs, 1:nkeep), &
init_overlaps => kp_init_overlaps(1:nkeep))
kp_final_hamil = 0.0_dp
do n = 1, nkeep
do m = 1, nkeep
do i = 1, n
do j = 1, m
kp_final_hamil(n, m) = kp_final_hamil(n, m) + S(i, n) * kp_hamil_mean(i, j) * S(j, m)
end do
end do
end do
end do
call dsyev('V', 'U', nkeep, kp_final_hamil, nkeep, kp_hamil_eigv, work, lwork, info)
init_overlaps = kp_final_hamil(:, 1) * sqrt(kp_overlap_mean(1, 1))
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)') kp_hamil_eigv(i), init_overlaps(i)
end do
end associate
deallocate(kp_final_hamil)
deallocate(kp_hamil_eigv)
end do
deallocate(work)
close(temp_unit)
end subroutine find_and_output_gs_eigv