find_and_output_lowdin_eigv Subroutine

public subroutine find_and_output_lowdin_eigv(config_label, nvecs, overlap_matrix, hamil_matrix, npositive, all_evals, tOutput, spin_matrix, all_spin)


Type IntentOptional Attributes Name
integer, intent(in) :: config_label
integer, intent(in) :: nvecs
real(kind=dp), intent(in) :: overlap_matrix(:,:)
real(kind=dp), intent(in) :: hamil_matrix(:,:)
integer, intent(out) :: npositive
real(kind=dp), intent(out) :: all_evals(:,:)
logical, intent(in) :: tOutput
real(kind=dp), intent(in), optional :: spin_matrix(:,:)
real(kind=dp), intent(out), optional :: all_spin(:,:)


Source Code

    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 &
        end if

        do nkeep = 1, npositive

            allocate(kp_final_hamil(nkeep, 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

            if (allocated(rotated_spin)) deallocate(rotated_spin)

        end do


        if (tOutput) close(temp_unit)

    end subroutine find_and_output_lowdin_eigv