SolutionHandler_module.f90

Go to the documentation of this file.

! Copyright 2019
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-in (UKRmol+ suite).
!
!     UKRmol-in is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
!
!     UKRmol-in is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with  UKRmol-in (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
 
module solutionhandler_module
 
    use cdenprop_defs,             only: civect, maxnuc
    use const_gbl,                 only: stdout
    use consts_mpi_ci,             only: symtype_d2h, save_continuum_coeffs, save_l2_coeffs, max_neig
    use precisn,                   only: wp
    use mpi_gbl,                   only: master, mpi_xermsg, mpi_mod_isend, mpi_mod_recv
    use params,                    only: ccidata, c8stars, cblank, cpoly
    use scatci_routines,           only: movep, civio, movew
    use baseintegral_module,       only: baseintegral
    use diagonalizerresult_module, only: diagonalizerresult
    use options_module,            only: options
    use parallelization_module,    only: grid => process_grid
    use sweden_module,             only: swedenintegral
    use ukrmol_module,             only: ukrmolintegral
 
    implicit none
 
    type, extends(diagonalizerresult) :: solutionhandler
        integer               :: io_unit = 25
        integer               :: symmetry_type
        integer               :: num_eigenpairs
        real(wp), allocatable :: energy_shifts(:)
        real(wp)              :: core_energy
        real(wp)              :: energy_shift
        integer, allocatable  :: phase(:)
        integer               :: size_phase
        integer               :: vec_dimen
        integer               :: nciset
        integer               :: current_eigenvector = 0
        integer               :: print_all_eigs = 0
    contains
        procedure, public :: construct
        procedure, public :: write_header       => write_header_sol
        procedure, public :: export_header      => export_header_sol
        procedure, public :: export_eigenvalues
        procedure, public :: shift_eigenvalues
        procedure, public :: handle_eigenvalues => write_eigenvalues
        procedure, public :: handle_eigenvector => write_eigenvector
        procedure, public :: finalize_solutions => finalize_solutions_sol
        procedure, public :: destroy
    end type solutionhandler
 
contains
 
    subroutine construct (this, opts)
 
        class(solutionhandler)      :: this
        class(options),  intent(in) :: opts
 
        integer :: n
 
        this % symmetry_type   = opts % sym_group_flag
        this % io_unit         = opts % ci_output_unit
        this % continuum_dimen = opts % seq_num_last_continuum_csf
        this % vec_dimen       = opts % contracted_mat_size
        this % num_eigenpairs  = opts % num_eigenpairs
        this % vector_storage  = opts % vector_storage_method
        this % print_all_eigs  = opts % print_flags(2)
 
        if (.not.(allocated(this % energy_shifts))) then
            allocate(this % energy_shifts(this % num_eigenpairs))
            n = min(this % num_eigenpairs, size(opts % energy_shifts))
            this % energy_shifts = 0.0_wp
            this % energy_shifts(1:n) = opts % energy_shifts(1:n)
        end if
 
    end subroutine construct
 
 
    subroutine destroy (this)
 
        class(solutionhandler), intent(inout) :: this
 
        if (allocated(this % energy_shifts)) then
            deallocate (this % energy_shifts)
        end if
 
        if (allocated(this % phase)) then
            deallocate (this % phase)
        end if
 
        call this % ci_vec % final_CV
 
    end subroutine destroy
 
 
    subroutine write_header_sol (this, option, integrals)
 
        class(solutionhandler)          :: this
        class(options),      intent(in) :: option
        class(baseintegral), intent(in) :: integrals
        integer                         :: Nset = 1
        integer                         :: print_header
        integer                         :: nnuc, nrec, i, buf(1), n = 1, tag = 0
        integer                         :: ierror, nth, nhd(10)
 
        !first lets quickly grab the core_energy while we're here
        this % core_energy = integrals % get_core_energy()
 
        if (grid % grank /= master) return
        this % nciset = option % output_ci_set_number
        if (option % output_ci_set_number < 0) return
 
        ! If output unit dataset for these solutions follows after another dataset in the same unit, we need to wait
        ! for a signal from master of the MPI group responsible for writing of that dataset first, so that we do not
        ! access the unit concurrently.
        if (option % preceding_writer >= 0) then
            call mpi_mod_recv(option % preceding_writer, tag, buf, n)
        end if
 
        this % size_phase = min(size(option % phase_index), this % vec_dimen)
 
        if (.not.(allocated(this % phase))) then
            allocate(this % phase(this % vec_dimen))
            this % phase = 0
            this % phase(1:this % size_phase) = option % phase_index(1:this % size_phase)
        end if
 
        nset = option % output_ci_set_number
        nth = nset
 
        print_header = option % print_dataset_heading
 
        if (this % symmetry_type >= symtype_d2h) then
 
            ! Move to the correct position in the file. The subroutine movep can, in theory, change the value of "nth". However,
            ! OptionsSet should have made sure that the values of output_ci_set_number are such that this will never happen.
            ! (I.e. that we are writing always the next unit in line.) Otherwise we would risk locking of the MPI groups and/or
            ! overwriting some MPI group results by other MPI groups.
            call movep(this % io_unit, nth, ierror, print_header, stdout)
 
            if (ierror /= 0) then
                call mpi_xermsg('SolutionHandler_module', 'write_header_sol', &
                                ' ERROR POSITIONING FOR OUTPUT OF CI COEFFICIENTS', 1, 1)
            end if
 
            nset = nth
            nnuc = integrals % get_num_nuclei()
            nrec = nnuc + this % num_eigenpairs + 1
 
            write (this % io_unit) c8stars, cblank, cblank, ccidata
            write (this % io_unit) nset, nrec, option % diag_name(1:120), nnuc,  option % contracted_mat_size, &
                                   this % num_eigenpairs, option % lambda, option % spin, option % spin_z, &
                                   option % num_electrons, this % core_energy, option % num_target_sym, &
                                   (option % lambda_continuum_orbitals_target(i), &
                                    option % num_continuum_orbitals_target(i), i = 1, option % num_target_sym)
 
            call integrals % write_geometries(this % io_unit)
 
        else
 
            ! Move to the correct position in the file
            call movew(this % io_unit, nth, ierror, print_header, stdout)
 
            if (ierror /= 0) then
                call mpi_xermsg('SolutionHandler_module', 'write_header_sol', &
                                ' ERROR POSITIONING FOR OUTPUT OF CI COEFFICIENTS', 2, 1)
            end if
 
            nhd( :) = 0
            nhd( 2) = option % contracted_mat_size
            nhd( 3) = this % num_eigenpairs
            nhd( 4) = option % num_syms
            nhd( 8) = option % contracted_mat_size
            nhd( 9) = integrals % get_num_nuclei()
            nhd(10) = merge(1, 0, this % num_eigenpairs < option % contracted_mat_size)
 
            ! We are using alchemy integrals
            write (this % io_unit) nth, nhd, option % diag_name(1:120), integrals % NHE, integrals % DTNUC
 
        end if
 
        write (stdout, "(' CI data will be stored as set number',I3)") nset
 
    end subroutine write_header_sol
 
 
    subroutine export_header_sol (this, option, integrals)
 
        class(solutionhandler)           :: this
        class(options),      intent(in)  :: option
        class(baseintegral), intent(in)  :: integrals
 
        integer :: err, i
 
        !first lets quickly grab the core_energy while we're here
        this % core_energy = integrals % get_core_energy()
        this % size_phase = min(size(option % phase_index), this % vec_dimen)
 
        if (.not.(allocated(this % phase))) then
            allocate(this % phase(this % vec_dimen), stat = err)
            if (err /= 0) then
                call mpi_xermsg('SolutionHandler_module', 'export_header_sol', &
                                'Error allocating internal phase array', 1, 1)
            end if
            this % phase = 0
            this % phase(1:this % size_phase) = option % phase_index(1:this % size_phase)
        end if
 
        if (this % symmetry_type >= symtype_d2h) then
            this % ci_vec % Nset  = option % output_ci_set_number
            this % ci_vec % nrec  = integrals % get_num_nuclei() + this % num_eigenpairs + 1
            this % ci_vec % name  = option % diag_name
            this % ci_vec % nnuc  = integrals % get_num_nuclei()
            this % ci_vec % nocsf = option % contracted_mat_size
            this % ci_vec % nstat = this % num_eigenpairs
            this % ci_vec % mgvn  = option % lambda
            this % ci_vec % s     = option % spin
            this % ci_vec % sz    = option % spin_z
            this % ci_vec % nelt  = option % num_electrons
            this % ci_vec % e0    = this % core_energy
 
            if (this % ci_vec % nnuc > maxnuc) then
                call mpi_xermsg('SolutionHandler_module', 'export_header_sol', &
                                'Increase the value of maxnuc in cdenprop_defs and recompile', 2, 1)
            end if
 
            select type (integrals)
                type is (ukrmolintegral)
                    this % ci_vec % cname (1:this % ci_vec % nnuc) = integrals % cname (1:this % ci_vec % nnuc)
                    this % ci_vec % charge(1:this % ci_vec % nnuc) = integrals % charge(1:this % ci_vec % nnuc)
                    this % ci_vec % xnuc  (1:this % ci_vec % nnuc) = integrals % xnuc  (1:this % ci_vec % nnuc)
                    this % ci_vec % ynuc  (1:this % ci_vec % nnuc) = integrals % ynuc  (1:this % ci_vec % nnuc)
                    this % ci_vec % znuc  (1:this % ci_vec % nnuc) = integrals % znuc  (1:this % ci_vec % nnuc)
                type is (swedenintegral)
                    this % ci_vec % cname (1:this % ci_vec % nnuc) = integrals % cname (1:this % ci_vec % nnuc)
                    this % ci_vec % charge(1:this % ci_vec % nnuc) = integrals % charge(1:this % ci_vec % nnuc)
                    this % ci_vec % xnuc  (1:this % ci_vec % nnuc) = integrals % xnuc  (1:this % ci_vec % nnuc)
                    this % ci_vec % ynuc  (1:this % ci_vec % nnuc) = integrals % ynuc  (1:this % ci_vec % nnuc)
                    this % ci_vec % znuc  (1:this % ci_vec % nnuc) = integrals % znuc  (1:this % ci_vec % nnuc)
                class default
                    call mpi_xermsg('SolutionHandler_module', 'export_header_sol', &
                                    'Geometry export implemented only for UKRMOLIntegral, SWEDENIntegral', 3, 1)
            end select
 
            write (stdout, '(/,"CI VECTOR HEADER DATA:")')
            write (stdout, '("----------------------")')
            write (stdout, *) this % ci_vec % nset,  &
                              this % ci_vec % nrec,  &
                              this % ci_vec % NAME,  &
                              this % ci_vec % nnuc,  &
                              this % ci_vec % nocsf, &
                              this % ci_vec % nstat, &
                              this % ci_vec % mgvn,  &
                              this % ci_vec % s,     &
                              this % ci_vec % sz,    &
                              this % ci_vec % nelt,  &
                              this % ci_vec % e0
            do i = 1, this % ci_vec % nnuc
                write (stdout, *) this % ci_vec % cname(i), &
                                  this % ci_vec % xnuc(i),  &
                                  this % ci_vec % ynuc(i),  &
                                  this % ci_vec % znuc(i),  &
                                  this % ci_vec % charge(i)
            end do
            write (stdout, '("----------------------")')
        else
            call mpi_xermsg('SolutionHandler_module', 'export_header_sol', &
                            'Output of ALCHEMY header not implemented yet', 4, 1)
        end if
 
        write (stdout, '(" CI data will be stored in memory for use in subsequent subroutines")')
 
    end subroutine export_header_sol
 
 
    subroutine shift_eigenvalues (this, output_eigenvalues, num_eigenpairs)
 
        class(solutionhandler)  :: this
        integer,  intent(in)    :: num_eigenpairs
        real(wp), intent(inout) :: output_eigenvalues(num_eigenpairs)
        integer :: ido, jdo, neig
 
        if (grid % grank /= master) return
        if (this % nciset < 0) return
 
        output_eigenvalues = output_eigenvalues + this % energy_shifts
 
        neig = this % num_eigenpairs
        if (this % print_all_eigs == 0) then
            neig = min(this % num_eigenpairs, max_neig)
        end if
 
        write (stdout, '("ENERGY SHIFTS TO BE APPLIED")')
        do ido = 1, neig, 5
            write (stdout, "(5F20.10)") (this % energy_shifts(jdo), jdo = ido, min(ido + 4, neig))
        end do
 
    end subroutine shift_eigenvalues
 
 
    subroutine write_eigenvalues (this, eigenvalues, diagonals, num_eigenpairs, vec_dimen)
 
        class(solutionhandler) :: this
        integer,  intent(in)   :: num_eigenpairs, vec_dimen
        real(wp), intent(in)   :: eigenvalues(num_eigenpairs), diagonals(vec_dimen)
        real(wp)               :: output_eigenvalues(num_eigenpairs)
        integer :: ido, jdo, neig
 
        if (grid % grank /= master) return
        if (this % nciset < 0) return
 
        output_eigenvalues = eigenvalues
 
        neig = this % num_eigenpairs
        if (this % print_all_eigs == 0) then
            neig = min(this % num_eigenpairs, max_neig)
        end if
 
        if (sum(abs(this % energy_shifts)) .gt. 0_wp) then
            write (stdout, '("ORIGINAL EIGEN-ENERGIES")')
            do ido = 1, neig, 5
                write (stdout, "(16X,5F20.10)") (output_eigenvalues(jdo) + this % core_energy, jdo = ido, min(ido + 4, neig))
            end do
            call this % shift_eigenvalues(output_eigenvalues, num_eigenpairs)
        end if
 
        write (stdout, '("EIGEN-ENERGIES")')
        do ido = 1, neig, 5
            write (stdout, "(16X,5F20.10)") (output_eigenvalues(jdo) + this % core_energy, jdo = ido, min(ido + 4, neig))
        end do
        write (this % io_unit) this % phase, output_eigenvalues, diagonals
 
    end subroutine write_eigenvalues
 
 
    subroutine export_eigenvalues (this, eigenvalues, diagonals, num_eigenpairs, vec_dimen, ei, iphz)
 
        class(solutionhandler) :: this
        integer,  intent(in)   :: num_eigenpairs, vec_dimen
        real(wp), intent(in)   :: eigenvalues(num_eigenpairs), diagonals(vec_dimen)
        real(wp), allocatable  :: ei(:)
        real(wp)               :: output_eigenvalues(num_eigenpairs)
        integer,  allocatable  :: iphz(:)
        integer                :: ido, jdo, err
 
        output_eigenvalues = eigenvalues
 
        if (sum(abs(this % energy_shifts)) .gt. 0_wp) then
            call this % shift_eigenvalues(output_eigenvalues, num_eigenpairs)
        end if
 
        if (allocated(ei)) deallocate(ei)
        allocate(ei, source = output_eigenvalues, stat = err)
        if (err /= 0) then
            call mpi_xermsg('SolutionHandler_module', 'export_eigenvalues', 'Error allocating ei', size(output_eigenvalues), 1)
        end if
 
        if (allocated(iphz)) deallocate(iphz)
        allocate(iphz, source = this % phase, stat = err)
        if (err /= 0) then
            call mpi_xermsg('SolutionHandler_module', 'export_eigenvalues', 'Error allocating iphz', size(this % phase), 1)
        end if
 
    end subroutine export_eigenvalues
 
 
    subroutine write_eigenvector (this, eigenvector, vec_dimen)
 
        class(solutionhandler) :: this
        integer,  intent(in)   :: vec_dimen
        real(wp), intent(in)   :: eigenvector(vec_dimen)
        integer                :: start_dimen, end_dimen, ido, j
 
        if (grid % grank /= master) return
        if (this % nciset < 0) return
 
        this % current_eigenvector = this % current_eigenvector + 1
 
        start_dimen = 1
        end_dimen = this % vec_dimen
 
        if (iand(this % vector_storage, save_l2_coeffs) == 0) then
            end_dimen = this % continuum_dimen
        end if
        if (iand(this % vector_storage, save_continuum_coeffs) == 0) then
            start_dimen = this % continuum_dimen + 1
        end if
 
        write (this % io_unit) this % current_eigenvector, (eigenvector(j), j = start_dimen, end_dimen)
 
    end subroutine write_eigenvector
 
 
    subroutine finalize_solutions_sol (this, option)
 
        class(solutionhandler)          :: this
        class(options),      intent(in) :: option
 
        integer :: buf(1), tag = 0, n = 1
 
        flush(this % io_unit)
        close(this % io_unit)
 
        if (option % following_writer >= 0) then
            call mpi_mod_isend(option % following_writer, buf, tag, n)
        end if
 
    end subroutine finalize_solutions_sol
 
end module solutionhandler_module