/scratch.ssd/codes/ukrmol-in-git/source/mpi-ci-diag/mpi_scatci.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/>.
 
!> \brief   MPI-SCATCI main program
!> \authors A Al-Refaie, J Benda
!> \date    2017 - 2022
!>
!> MPI-SCATCI diagonalizes Hamiltonians for one or more irreducible representations. The input for each
!> is read from the single input file as the pair of standard SCATCI namelists &input and &cinorn.
!> If there are fewer MPI processes than Hamiltonians to diagonalize, each process is given several
!> Hamiltonians to diagonalize sequentially. If there are more processes than Hamiltonians, some or all
!> Hamiltonians will be diagonalized in a distributed regime. All processes that participate on diagonalization
!> of the same matrix are part of a "diagonalization group" which has separate BLACS context and their own
!> MPI communicator.
!>
!> The eigenvectors are kept in the memory and, if desired, are then used to calculate transition dipole
!> moments and further derived data.
!>
!> \note 16/01/2019 - Jakub Benda: Unifom coding style and expanded documentation.
!> \note 25/01/2019 - Jakub Benda: Extension to support multiple subsequent diagonalizations.
!> \note 30/01/2019 - Jakub Benda: Evaluation of properties using libcdenprop, including cross-symmetry ones.
!> \note 22/05/2022 - Jakub Benda: Reimplemented reading of eigensolutions from disk.
!>
program mpi_scatci
 
    use consts_mpi_ci,                   only: no_diagonalization, main_hamiltonian, target_hamiltonian, no_ci_target, &
                                               read_from_file
    use const_gbl,                       only: stdout
    use global_utils,                    only: print_ukrmol_header
    use precisn,                         only: wp, longint
    use mpi_gbl,                         only: mpi_mod_start, mpi_mod_finalize, master, myrank, mpi_mod_print_info
    use options_module,                  only: optionsset, ham_pars_io_wrapper
    use orbital_module,                  only: orbitaltable
    use csf_module,                      only: csfobject, csforbital, csfmanager
    use ci_hamiltonian_module,           only: target_ci_hamiltonian
    use basematrix_module,               only: basematrix
    use baseintegral_module,             only: baseintegral
    use diagonalizer_module,             only: basediagonalizer
    use memorymanager_module,            only: master_memory
    use parallelization_module,          only: process_grid
    use postprocessing_module,           only: postprocess
    use uncontracted_hamiltonian_module, only: uncontracted_hamiltonian
    use contracted_hamiltonian_module,   only: contracted_hamiltonian
    use timing_module,                   only: master_timer
    use target_rmat_ci_module,           only: target_rmat_ci, read_ci_mat
    use dispatcher_module,               only: dispatchintegral, dispatchmatrixanddiagonalizer, initialize_libraries
    use writermatrix_module,             only: writermatrix
    use solutionhandler_module,          only: solutionhandler
 
    implicit none
 
    type(optionsset)                        ::      scatci_input
    type(orbitaltable)                      ::      orbitals
    type(uncontracted_hamiltonian)          ::      enrgms
    type(contracted_hamiltonian)            ::      enrgmx
    type(csfmanager)                        ::      configuration_manager
    class(baseintegral),          pointer   ::      integrals           => null()
    class(target_ci_hamiltonian), pointer   ::      tgt_ci_hamiltonian  => null()
    class(basediagonalizer),      pointer   ::      diagonalizer        => null()
    class(basematrix),            pointer   ::      matrix_elements     => null()
    type(solutionhandler),    allocatable   ::      solutions(:)
    type(csfobject),          allocatable   ::      csfs(:)
    type(target_rmat_ci),     allocatable   ::      ci_rmat(:)
    real(wp),                 allocatable   ::      test_eig(:), test_vecs(:,:)
    integer                                 ::      num_mat_elms, i, j
    real(wp)                                ::      usr_ef
    integer                                 ::      usr_ef_type
 
    logical, parameter :: sequential_diagonalizations = .false.
    logical, parameter :: master_writes_to_stdout = .true.
    logical, parameter :: allow_shared_memory = .true.
 
    call mpi_mod_start(master_writes_to_stdout, allow_shared_memory)
    call master_timer % initialize()
    if (myrank == master) call print_ukrmol_header(stdout)
    call mpi_mod_print_info(stdout)
 
    ! Initialize all libraries used by MPI-SCATCI that need it (e.g. SLEPc)
    call initialize_libraries
 
    ! Read options (for all symmetries)
    call scatci_input % read
    allocate (solutions(size(scatci_input % opts)))
 
    ! Set up the MPI/BLACS groups
    call process_grid % setup(size(scatci_input % opts), sequential_diagonalizations)
    call process_grid % summarize
    if (.not. process_grid % sequential) call scatci_input % setup_write_order
 
    call master_timer % start_timer("Wall Time")
 
    ! Process all symmetries
    symmetry_loop: do i = 1, size(scatci_input % opts)
 
        if (.not. process_grid % is_my_group_work(i)) cycle
 
        call orbitals % initialize(scatci_input % opts(i) % tot_num_spin_orbitals,    &
                                   scatci_input % opts(i) % tot_num_orbitals,         &
                                   scatci_input % opts(i) % sym_group_flag,           &
                                   scatci_input % opts(i) % num_syms,                 &
                                   scatci_input % opts(i) % positron_flag)
        call orbitals % construct(scatci_input % opts(i) % num_orbitals_sym,                    &
                                  scatci_input % opts(i) % num_target_orbitals_sym_dinf_congen, &
                                  scatci_input % opts(i) % num_orbitals_sym,                    &
                                  scatci_input % opts(i) % num_electronic_orbitals_sym,         &
                                  scatci_input % opts(i) % num_target_orbitals_sym_congen)
 
        ! Compute the electron number for each reference orbital
        call orbitals % compute_electron_index(scatci_input % opts(i) % num_electrons, &
                                               scatci_input % opts(i) % reference_dtrs)
 
        ! Construct our csf manager
        call configuration_manager % initialize(scatci_input % opts(i), orbitals)
 
        ! Construct our CSFS
        call master_timer % start_timer("Construct CSFs")
        call configuration_manager % create_csfs(csfs,                      &
                scatci_input % opts(i) % orbital_sequence_number,           &
                scatci_input % opts(i) % num_ci_target_sym,                 &
                scatci_input % opts(i) % lambda_continuum_orbitals_target)
        call master_timer % stop_timer("Construct CSFs")
        call master_memory % print_memory_report
 
        ! Read molecular integrals (assuming all symmetries use the same!)
        if (.not. associated(integrals)) then
            call dispatchintegral(scatci_input % opts(i) % sym_group_flag, &
                                  scatci_input % opts(i) % use_UKRMOL_integrals, &
                                  integrals)
            call integrals % initialize(scatci_input % opts(i), &
                                        orbitals % total_num_orbitals, &
                                        orbitals % orbital_map)
            call integrals % load_integrals(scatci_input % opts(i) % integral_unit)
            call master_memory % print_memory_report
        end if
 
        ! if only reading already calculated eigensolutions (for processing in CDENPROP or outer interface)
        if (scatci_input % opts(i) % diagonalization_flag == read_from_file) then
            call solutions(i) % construct(scatci_input % opts(i))
            call solutions(i) % export_header(scatci_input % opts(i), integrals)
            call solutions(i) % read_from_file(scatci_input % opts(i))
            deallocate(csfs)
            call configuration_manager % finalize
            cycle
        end if
 
        write (stdout, *) '  Load CI target and construct matrix elements'
 
        ! Before we do anything else let us get our ECP if we need them
       !call DispatchECP(SCATCI_input % opts(i) % ecp_type,             &
       !                 SCATCI_input % opts(i) % ecp_filename,         &
       !                 SCATCI_input % opts(i) % all_ecp_defined,      &
       !                 SCATCI_input % opts(i) % num_target_sym,       &
       !                 SCATCI_input % opts(i) % num_target_state_sym, &
       !                 SCATCI_input % opts(i) % target_spatial,       &
       !                 SCATCI_input % opts(i) % target_multiplicity, ecp)
 
        if (scatci_input % opts(i) % do_ci_contraction()) then
            usr_ef = scatci_input % opts(i) % enh_factor
            usr_ef_type = scatci_input % opts(i) % enh_factor_type
            scatci_input % opts(i) % enh_factor = 1.0
            scatci_input % opts(i) % enh_factor_type = 0
            ! allocate the CIRmats
            allocate(ci_rmat(scatci_input % opts(i) % num_target_sym))
            do j = 1, scatci_input % opts(i) % num_target_sym
                call ci_rmat(j) % initialize(j,                                                &
                                             scatci_input % opts(i) % num_target_state_sym(j), &
                                             scatci_input % opts(i) % num_ci_target_sym(j),    &
                                             scatci_input % opts(i) % ci_phase(j),             &
                                             integrals % get_core_energy())
            end do
 
            if (scatci_input % opts(i) % ci_target_switch > 0) then
                call read_ci_mat(scatci_input % opts(i), ci_rmat)
            else
                do j = 1, scatci_input % opts(i) % num_target_sym
                    allocate(target_ci_hamiltonian::tgt_ci_hamiltonian)
                    call dispatchmatrixanddiagonalizer(scatci_input % opts(i) % diagonalizer_choice,    &
                                                       scatci_input % opts(i) % force_serial,           &
                                                       scatci_input % opts(i) % num_ci_target_sym(j),   &
                                                       ci_rmat(j) % nstat,                              &
                                                       matrix_elements,                                 &
                                                       diagonalizer,                                    &
                                                       integrals,                                       &
                                                       scatci_input % opts(i) % hamiltonian_unit)
 
                    call matrix_elements % construct(target_hamiltonian + j)
                    call matrix_elements % set_options(scatci_input % opts(i))
 
                    call tgt_ci_hamiltonian % construct(scatci_input % opts(i), csfs, orbitals, integrals)
                    call tgt_ci_hamiltonian % initialize(j)
 
                    call master_timer % start_timer("Build CI Hamiltonian")
                    call tgt_ci_hamiltonian % build_hamiltonian(matrix_elements)
                    call master_timer % stop_timer("Build CI Hamiltonian")
                    call master_timer % report_timers
 
                    call master_timer % start_timer("Diagonalization")
                    call diagonalizer % diagonalize(matrix_elements, ci_rmat(j) % nstat, ci_rmat(j), .true., &
                                                    scatci_input % opts(i), integrals)
                    call master_timer % stop_timer("Diagonalization")
 
                    call ci_rmat(j) % print
                    call master_timer % report_timers
 
                    call matrix_elements % destroy
                    if (process_grid % grank == master) close (scatci_input % opts(i) % hamiltonian_unit)
                    deallocate(matrix_elements, diagonalizer)
                    deallocate(tgt_ci_hamiltonian)
                end do
            end if
            scatci_input % opts(i) % enh_factor = usr_ef
            scatci_input % opts(i) % enh_factor_type = usr_ef_type
        end if
 
        if (scatci_input % opts(i) % diagonalization_flag /= no_diagonalization) then
            call dispatchmatrixanddiagonalizer(scatci_input % opts(i) % diagonalizer_choice,  &
                                               scatci_input % opts(i) % force_serial,         &
                                               scatci_input % opts(i) % contracted_mat_size,  &
                                               scatci_input % opts(i) % num_eigenpairs,       &
                                               matrix_elements,                               &
                                               diagonalizer,                                  &
                                               integrals,                                     &
                                               scatci_input % opts(i) % hamiltonian_unit)
            call matrix_elements % construct(main_hamiltonian)
            call matrix_elements % set_options(scatci_input % opts(i))
        else
            call integrals % write_matrix_header(scatci_input % opts(i) % hamiltonian_unit, &
                                                 scatci_input % opts(i) % contracted_mat_size)
            allocate(writermatrix::matrix_elements)
            call matrix_elements % construct(main_hamiltonian)
            call matrix_elements % set_options(scatci_input % opts(i))
        end if
 
        call matrix_elements % exclude_row_column(scatci_input % opts(i) % exclude_rowcolumn)
 
        ! We are doing a scattering calculation
        if (scatci_input % opts(i) % do_ci_contraction()) then
            call enrgmx % construct(scatci_input % opts(i), csfs, orbitals, integrals)
            call enrgmx % initialize(ci_rmat)
            call master_timer % start_timer("C-Hamiltonian Build")
            call enrgmx % build_hamiltonian(matrix_elements)
            call master_timer % stop_timer("C-Hamiltonian Build")
            deallocate (ci_rmat)
        else
            call enrgms % construct(scatci_input % opts(i), csfs, orbitals, integrals)
            call master_timer % start_timer("Target-Hamiltonian Build")
            call enrgms % build_hamiltonian(matrix_elements)
            call master_timer % stop_timer("Target-Hamiltonian Build")
        end if
 
        ! Lets deallocate all the csfs as well
        deallocate(csfs)
 
        ! Clear the integrals, we dont need them anymore
        call configuration_manager % finalize
 
        call master_memory % print_memory_report
        write (stdout, "('Matrix N is ',i8)") matrix_elements % get_matrix_size()
        write (stdout, "('Num elements is ',i14)") matrix_elements % get_size()
        call master_timer % report_timers
 
        if (scatci_input % opts(i) % diagonalization_flag /= no_diagonalization) then
            scatci_input % opts(i) % num_eigenpairs = min(scatci_input % opts(i) % num_eigenpairs, &
                                                          matrix_elements % get_matrix_size())
            call solutions(i) % construct(scatci_input % opts(i))
            call master_timer % start_timer("Diagonalization")
            call diagonalizer % diagonalize(matrix_elements,                          &
                                            scatci_input % opts(i) % num_eigenpairs,  &
                                            solutions(i),                             &
                                            .false.,                                  &
                                            scatci_input % opts(i),                   &
                                            integrals)
            call master_timer % stop_timer("Diagonalization")
            call master_timer % report_timers
        else
            write (stdout, '(/,"Diagonalization not selected")')
            num_mat_elms = matrix_elements % get_size()
 
            call master_timer % start_timer("Matrix Save")
           !call matrix_elements % save(SCATCI_input % hamiltonian_unit, SCATCI_input % num_matrix_elements_per_rec)
            call master_timer % stop_timer("Matrix Save")
 
            if (process_grid % grank == master) call ham_pars_io_wrapper(scatci_input % opts(i), .true., num_mat_elms)
 
            write (stdout, '(/,"Parameters written to: ham_data")')
        end if
 
        if (process_grid % grank == master) close (scatci_input % opts(i) % hamiltonian_unit)
 
        call matrix_elements % destroy
 
        deallocate (diagonalizer)
        deallocate (matrix_elements)
 
    end do symmetry_loop
 
    ! use solutions in outer interface etc.
    call postprocess(scatci_input, solutions)
    call master_timer % stop_timer("Wall Time")
    call master_timer % report_timers
 
    ! release (possibly shared) memory occupied by the integral storage
    if (associated(integrals)) then
        call integrals % finalize
        nullify (integrals)
    end if
 
    ! cleanly exit MPI
    call mpi_mod_finalize
 
end program mpi_scatci