LCOV - code coverage report
Current view: top level - src - qs_mo_io.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:4dc10b3) Lines: 537 594 90.4 %
Date: 2024-11-21 06:45:46 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Definition and initialisation of the mo data type.
      10             : !> \par History
      11             : !>      - adapted to the new QS environment data structure (02.04.2002,MK)
      12             : !>      - set_mo_occupation added (17.04.02,MK)
      13             : !>      - correct_mo_eigenvalues added (18.04.02,MK)
      14             : !>      - calculate_density_matrix moved from qs_scf to here (22.04.02,MK)
      15             : !>      - mo_set_p_type added (23.04.02,MK)
      16             : !>      - PRIVATE attribute set for TYPE mo_set_type (23.04.02,MK)
      17             : !>      - started conversion to LSD (1.2003, Joost VandeVondele)
      18             : !>      - Split of from qs_mo_types (07.2014, JGH)
      19             : !> \author Matthias Krack (09.05.2001,MK)
      20             : ! **************************************************************************************************
      21             : MODULE qs_mo_io
      22             : 
      23             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      24             :                                               get_atomic_kind,&
      25             :                                               get_atomic_kind_set
      26             :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      27             :                                               gto_basis_set_type
      28             :    USE cp_dbcsr_api,                    ONLY: dbcsr_binary_write,&
      29             :                                               dbcsr_checksum,&
      30             :                                               dbcsr_create,&
      31             :                                               dbcsr_p_type,&
      32             :                                               dbcsr_release,&
      33             :                                               dbcsr_type
      34             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      35             :                                               copy_fm_to_dbcsr
      36             :    USE cp_files,                        ONLY: close_file,&
      37             :                                               open_file
      38             :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      39             :                                               cp_fm_get_submatrix,&
      40             :                                               cp_fm_set_all,&
      41             :                                               cp_fm_set_submatrix,&
      42             :                                               cp_fm_to_fm,&
      43             :                                               cp_fm_type,&
      44             :                                               cp_fm_write_unformatted
      45             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      46             :                                               cp_logger_get_default_unit_nr,&
      47             :                                               cp_logger_type,&
      48             :                                               cp_to_string
      49             :    USE cp_output_handling,              ONLY: cp_p_file,&
      50             :                                               cp_print_key_finished_output,&
      51             :                                               cp_print_key_generate_filename,&
      52             :                                               cp_print_key_should_output,&
      53             :                                               cp_print_key_unit_nr
      54             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      55             :                                               section_vals_type,&
      56             :                                               section_vals_val_get
      57             :    USE kahan_sum,                       ONLY: accurate_sum
      58             :    USE kinds,                           ONLY: default_path_length,&
      59             :                                               default_string_length,&
      60             :                                               dp
      61             :    USE message_passing,                 ONLY: mp_para_env_type
      62             :    USE orbital_pointers,                ONLY: indco,&
      63             :                                               nco,&
      64             :                                               nso
      65             :    USE orbital_symbols,                 ONLY: cgf_symbol,&
      66             :                                               sgf_symbol
      67             :    USE orbital_transformation_matrices, ONLY: orbtramat
      68             :    USE particle_types,                  ONLY: particle_type
      69             :    USE physcon,                         ONLY: evolt
      70             :    USE qs_density_matrices,             ONLY: calculate_density_matrix
      71             :    USE qs_dftb_types,                   ONLY: qs_dftb_atom_type
      72             :    USE qs_dftb_utils,                   ONLY: get_dftb_atom_param
      73             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      74             :                                               get_qs_kind_set,&
      75             :                                               qs_kind_type
      76             :    USE qs_mo_occupation,                ONLY: set_mo_occupation
      77             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      78             :                                               mo_set_type
      79             : #include "./base/base_uses.f90"
      80             : 
      81             :    IMPLICIT NONE
      82             : 
      83             :    PRIVATE
      84             : 
      85             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_io'
      86             : 
      87             :    PUBLIC :: wfn_restart_file_name, &
      88             :              write_rt_mos_to_restart, &
      89             :              read_rt_mos_from_restart, &
      90             :              write_dm_binary_restart, &
      91             :              write_mo_set_to_output_unit, &
      92             :              write_mo_set_to_restart, &
      93             :              read_mo_set_from_restart, &
      94             :              read_mos_restart_low, &
      95             :              write_mo_set_low
      96             : 
      97             : CONTAINS
      98             : 
      99             : ! **************************************************************************************************
     100             : !> \brief ...
     101             : !> \param mo_array ...
     102             : !> \param particle_set ...
     103             : !> \param dft_section ...
     104             : !> \param qs_kind_set ...
     105             : ! **************************************************************************************************
     106      143589 :    SUBROUTINE write_mo_set_to_restart(mo_array, particle_set, dft_section, qs_kind_set)
     107             : 
     108             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     109             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     110             :       TYPE(section_vals_type), POINTER                   :: dft_section
     111             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     112             : 
     113             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_mo_set_to_restart'
     114             :       CHARACTER(LEN=30), DIMENSION(2), PARAMETER :: &
     115             :          keys = (/"SCF%PRINT%RESTART_HISTORY", "SCF%PRINT%RESTART        "/)
     116             : 
     117             :       INTEGER                                            :: handle, ikey, ires, ispin
     118             :       TYPE(cp_logger_type), POINTER                      :: logger
     119             : 
     120      143589 :       CALL timeset(routineN, handle)
     121      143589 :       logger => cp_get_default_logger()
     122             : 
     123             :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     124      143589 :                                            dft_section, keys(1)), cp_p_file) .OR. &
     125             :           BTEST(cp_print_key_should_output(logger%iter_info, &
     126             :                                            dft_section, keys(2)), cp_p_file)) THEN
     127             : 
     128       15875 :          IF (mo_array(1)%use_mo_coeff_b) THEN
     129             :             ! we are using the dbcsr mo_coeff
     130             :             ! we copy it to the fm for anycase
     131       11316 :             DO ispin = 1, SIZE(mo_array)
     132        6135 :                IF (.NOT. ASSOCIATED(mo_array(ispin)%mo_coeff_b)) THEN
     133           0 :                   CPASSERT(.FALSE.)
     134             :                END IF
     135             :                CALL copy_dbcsr_to_fm(mo_array(ispin)%mo_coeff_b, &
     136       11316 :                                      mo_array(ispin)%mo_coeff) !fm->dbcsr
     137             :             END DO
     138             :          END IF
     139             : 
     140       47625 :          DO ikey = 1, SIZE(keys)
     141             : 
     142       31750 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     143      143589 :                                                  dft_section, keys(ikey)), cp_p_file)) THEN
     144             : 
     145             :                ires = cp_print_key_unit_nr(logger, dft_section, keys(ikey), &
     146             :                                            extension=".wfn", file_status="REPLACE", file_action="WRITE", &
     147       15887 :                                            do_backup=.TRUE., file_form="UNFORMATTED")
     148             : 
     149             :                CALL write_mo_set_low(mo_array, particle_set=particle_set, &
     150       15887 :                                      qs_kind_set=qs_kind_set, ires=ires)
     151             : 
     152       15887 :                CALL cp_print_key_finished_output(ires, logger, dft_section, TRIM(keys(ikey)))
     153             :             END IF
     154             :          END DO
     155             :       END IF
     156             : 
     157      143589 :       CALL timestop(handle)
     158             : 
     159      143589 :    END SUBROUTINE write_mo_set_to_restart
     160             : 
     161             : ! **************************************************************************************************
     162             : !> \brief calculates density matrix from mo set and writes the density matrix
     163             : !>        into a binary restart file
     164             : !> \param mo_array mos
     165             : !> \param dft_section dft input section
     166             : !> \param tmpl_matrix template dbcsr matrix
     167             : !> \author Mohammad Hossein Bani-Hashemian
     168             : ! **************************************************************************************************
     169        9487 :    SUBROUTINE write_dm_binary_restart(mo_array, dft_section, tmpl_matrix)
     170             : 
     171             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     172             :       TYPE(section_vals_type), POINTER                   :: dft_section
     173             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: tmpl_matrix
     174             : 
     175             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_dm_binary_restart'
     176             : 
     177             :       CHARACTER(LEN=default_path_length)                 :: file_name, project_name
     178             :       INTEGER                                            :: handle, ispin, unit_nr
     179             :       LOGICAL                                            :: do_dm_restart
     180             :       REAL(KIND=dp)                                      :: cs_pos
     181             :       TYPE(cp_logger_type), POINTER                      :: logger
     182             :       TYPE(dbcsr_type), POINTER                          :: matrix_p_tmp
     183             : 
     184        9487 :       CALL timeset(routineN, handle)
     185        9487 :       logger => cp_get_default_logger()
     186        9487 :       IF (logger%para_env%is_source()) THEN
     187        4871 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     188             :       ELSE
     189             :          unit_nr = -1
     190             :       END IF
     191             : 
     192        9487 :       project_name = logger%iter_info%project_name
     193        9487 :       CALL section_vals_val_get(dft_section, "SCF%PRINT%DM_RESTART_WRITE", l_val=do_dm_restart)
     194        9487 :       NULLIFY (matrix_p_tmp)
     195             : 
     196        9487 :       IF (do_dm_restart) THEN
     197           0 :          ALLOCATE (matrix_p_tmp)
     198           0 :          DO ispin = 1, SIZE(mo_array)
     199           0 :             CALL dbcsr_create(matrix_p_tmp, template=tmpl_matrix(ispin)%matrix, name="DM RESTART")
     200             : 
     201           0 :             IF (.NOT. ASSOCIATED(mo_array(ispin)%mo_coeff_b)) CPABORT("mo_coeff_b NOT ASSOCIATED")
     202             : 
     203           0 :             CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, mo_array(ispin)%mo_coeff_b)
     204             :             CALL calculate_density_matrix(mo_array(ispin), matrix_p_tmp, &
     205           0 :                                           use_dbcsr=.TRUE., retain_sparsity=.FALSE.)
     206             : 
     207           0 :             WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_SCF_DM_SPIN_", ispin, "_RESTART.dm"
     208           0 :             cs_pos = dbcsr_checksum(matrix_p_tmp, pos=.TRUE.)
     209           0 :             IF (unit_nr > 0) THEN
     210           0 :                WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
     211             :             END IF
     212           0 :             CALL dbcsr_binary_write(matrix_p_tmp, file_name)
     213             : 
     214           0 :             CALL dbcsr_release(matrix_p_tmp)
     215             :          END DO
     216           0 :          DEALLOCATE (matrix_p_tmp)
     217             :       END IF
     218             : 
     219        9487 :       CALL timestop(handle)
     220             : 
     221        9487 :    END SUBROUTINE write_dm_binary_restart
     222             : 
     223             : ! **************************************************************************************************
     224             : !> \brief ...
     225             : !> \param mo_array ...
     226             : !> \param rt_mos ...
     227             : !> \param particle_set ...
     228             : !> \param dft_section ...
     229             : !> \param qs_kind_set ...
     230             : ! **************************************************************************************************
     231         432 :    SUBROUTINE write_rt_mos_to_restart(mo_array, rt_mos, particle_set, dft_section, &
     232             :                                       qs_kind_set)
     233             : 
     234             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     235             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: rt_mos
     236             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     237             :       TYPE(section_vals_type), POINTER                   :: dft_section
     238             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     239             : 
     240             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_mos_to_restart'
     241             :       CHARACTER(LEN=43), DIMENSION(2), PARAMETER :: keys = (/ &
     242             :          "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY", &
     243             :          "REAL_TIME_PROPAGATION%PRINT%RESTART        "/)
     244             : 
     245             :       INTEGER                                            :: handle, ikey, ires
     246             :       TYPE(cp_logger_type), POINTER                      :: logger
     247             : 
     248         432 :       CALL timeset(routineN, handle)
     249         432 :       logger => cp_get_default_logger()
     250             : 
     251             :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     252         432 :                                            dft_section, keys(1)), cp_p_file) .OR. &
     253             :           BTEST(cp_print_key_should_output(logger%iter_info, &
     254             :                                            dft_section, keys(2)), cp_p_file)) THEN
     255             : 
     256         324 :          DO ikey = 1, SIZE(keys)
     257             : 
     258         216 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     259         432 :                                                  dft_section, keys(ikey)), cp_p_file)) THEN
     260             : 
     261             :                ires = cp_print_key_unit_nr(logger, dft_section, keys(ikey), &
     262             :                                            extension=".rtpwfn", file_status="REPLACE", file_action="WRITE", &
     263         108 :                                            do_backup=.TRUE., file_form="UNFORMATTED")
     264             : 
     265             :                CALL write_mo_set_low(mo_array, rt_mos=rt_mos, qs_kind_set=qs_kind_set, &
     266         108 :                                      particle_set=particle_set, ires=ires)
     267         108 :                CALL cp_print_key_finished_output(ires, logger, dft_section, TRIM(keys(ikey)))
     268             :             END IF
     269             :          END DO
     270             :       END IF
     271             : 
     272         432 :       CALL timestop(handle)
     273             : 
     274         432 :    END SUBROUTINE write_rt_mos_to_restart
     275             : 
     276             : ! **************************************************************************************************
     277             : !> \brief ...
     278             : !> \param mo_array ...
     279             : !> \param qs_kind_set ...
     280             : !> \param particle_set ...
     281             : !> \param ires ...
     282             : !> \param rt_mos ...
     283             : ! **************************************************************************************************
     284       16003 :    SUBROUTINE write_mo_set_low(mo_array, qs_kind_set, particle_set, ires, rt_mos)
     285             : 
     286             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     287             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     288             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     289             :       INTEGER                                            :: ires
     290             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
     291             :          OPTIONAL                                        :: rt_mos
     292             : 
     293             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'write_mo_set_low'
     294             : 
     295             :       INTEGER                                            :: handle, iatom, ikind, imat, iset, &
     296             :                                                             ishell, ispin, lmax, lshell, &
     297             :                                                             max_block, nao, natom, nmo, nset, &
     298             :                                                             nset_max, nshell_max, nspin
     299       16003 :       INTEGER, DIMENSION(:), POINTER                     :: nset_info, nshell
     300       16003 :       INTEGER, DIMENSION(:, :), POINTER                  :: l, nshell_info
     301       16003 :       INTEGER, DIMENSION(:, :, :), POINTER               :: nso_info
     302             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     303             :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
     304             : 
     305       16003 :       CALL timeset(routineN, handle)
     306       16003 :       nspin = SIZE(mo_array)
     307       16003 :       nao = mo_array(1)%nao
     308             : 
     309       16003 :       IF (ires > 0) THEN
     310             :          !     *** create some info about the basis set first ***
     311        8173 :          natom = SIZE(particle_set, 1)
     312        8173 :          nset_max = 0
     313        8173 :          nshell_max = 0
     314             : 
     315       52615 :          DO iatom = 1, natom
     316       44442 :             NULLIFY (orb_basis_set, dftb_parameter)
     317       44442 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     318             :             CALL get_qs_kind(qs_kind_set(ikind), &
     319       44442 :                              basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
     320       97057 :             IF (ASSOCIATED(orb_basis_set)) THEN
     321             :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     322             :                                       nset=nset, &
     323             :                                       nshell=nshell, &
     324       37008 :                                       l=l)
     325       37008 :                nset_max = MAX(nset_max, nset)
     326      104396 :                DO iset = 1, nset
     327      104396 :                   nshell_max = MAX(nshell_max, nshell(iset))
     328             :                END DO
     329        7434 :             ELSEIF (ASSOCIATED(dftb_parameter)) THEN
     330        7433 :                CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
     331        7433 :                nset_max = MAX(nset_max, 1)
     332        7433 :                nshell_max = MAX(nshell_max, lmax + 1)
     333             :             ELSE
     334             :                ! We assume here an atom without a basis set
     335             :                ! CPABORT("Unknown basis type. ")
     336             :             END IF
     337             :          END DO
     338             : 
     339       40865 :          ALLOCATE (nso_info(nshell_max, nset_max, natom))
     340      304713 :          nso_info(:, :, :) = 0
     341             : 
     342       32692 :          ALLOCATE (nshell_info(nset_max, natom))
     343      136333 :          nshell_info(:, :) = 0
     344             : 
     345       24519 :          ALLOCATE (nset_info(natom))
     346       52615 :          nset_info(:) = 0
     347             : 
     348       52615 :          DO iatom = 1, natom
     349       44442 :             NULLIFY (orb_basis_set, dftb_parameter)
     350       44442 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     351             :             CALL get_qs_kind(qs_kind_set(ikind), &
     352       44442 :                              basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
     353       97057 :             IF (ASSOCIATED(orb_basis_set)) THEN
     354             :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     355             :                                       nset=nset, &
     356             :                                       nshell=nshell, &
     357       37008 :                                       l=l)
     358       37008 :                nset_info(iatom) = nset
     359      104396 :                DO iset = 1, nset
     360       67388 :                   nshell_info(iset, iatom) = nshell(iset)
     361      200720 :                   DO ishell = 1, nshell(iset)
     362       96324 :                      lshell = l(ishell, iset)
     363      163712 :                      nso_info(ishell, iset, iatom) = nso(lshell)
     364             :                   END DO
     365             :                END DO
     366        7434 :             ELSEIF (ASSOCIATED(dftb_parameter)) THEN
     367        7433 :                CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
     368        7433 :                nset_info(iatom) = 1
     369        7433 :                nshell_info(1, iatom) = lmax + 1
     370       17412 :                DO ishell = 1, lmax + 1
     371        9979 :                   lshell = ishell - 1
     372       17412 :                   nso_info(ishell, 1, iatom) = nso(lshell)
     373             :                END DO
     374             :             ELSE
     375             :                ! We assume here an atom without a basis set
     376             :                ! CPABORT("Unknown basis type. ")
     377             :             END IF
     378             :          END DO
     379             : 
     380        8173 :          WRITE (ires) natom, nspin, nao, nset_max, nshell_max
     381       52615 :          WRITE (ires) nset_info
     382      136333 :          WRITE (ires) nshell_info
     383      304713 :          WRITE (ires) nso_info
     384             : 
     385        8173 :          DEALLOCATE (nset_info)
     386             : 
     387        8173 :          DEALLOCATE (nshell_info)
     388             : 
     389        8173 :          DEALLOCATE (nso_info)
     390             :       END IF
     391             : 
     392             :       ! use the scalapack block size as a default for buffering columns
     393       16003 :       CALL cp_fm_get_info(mo_array(1)%mo_coeff, ncol_block=max_block)
     394       34822 :       DO ispin = 1, nspin
     395       18819 :          nmo = mo_array(ispin)%nmo
     396       18819 :          IF ((ires > 0) .AND. (nmo > 0)) THEN
     397        9476 :             WRITE (ires) nmo, &
     398        9476 :                mo_array(ispin)%homo, &
     399        9476 :                mo_array(ispin)%lfomo, &
     400       18952 :                mo_array(ispin)%nelectron
     401       90275 :             WRITE (ires) mo_array(ispin)%eigenvalues(1:nmo), &
     402       99751 :                mo_array(ispin)%occupation_numbers(1:nmo)
     403             :          END IF
     404       34822 :          IF (PRESENT(rt_mos)) THEN
     405         414 :             DO imat = 2*ispin - 1, 2*ispin
     406         414 :                CALL cp_fm_write_unformatted(rt_mos(imat), ires)
     407             :             END DO
     408             :          ELSE
     409       18681 :             CALL cp_fm_write_unformatted(mo_array(ispin)%mo_coeff, ires)
     410             :          END IF
     411             :       END DO
     412             : 
     413       16003 :       CALL timestop(handle)
     414             : 
     415       16003 :    END SUBROUTINE write_mo_set_low
     416             : 
     417             : ! **************************************************************************************************
     418             : !> \brief ...
     419             : !> \param filename ...
     420             : !> \param exist ...
     421             : !> \param section ...
     422             : !> \param logger ...
     423             : !> \param kp ...
     424             : !> \param xas ...
     425             : !> \param rtp ...
     426             : ! **************************************************************************************************
     427        1184 :    SUBROUTINE wfn_restart_file_name(filename, exist, section, logger, kp, xas, rtp)
     428             :       CHARACTER(LEN=default_path_length), INTENT(OUT)    :: filename
     429             :       LOGICAL, INTENT(OUT)                               :: exist
     430             :       TYPE(section_vals_type), POINTER                   :: section
     431             :       TYPE(cp_logger_type), POINTER                      :: logger
     432             :       LOGICAL, INTENT(IN), OPTIONAL                      :: kp, xas, rtp
     433             : 
     434             :       INTEGER                                            :: n_rep_val
     435             :       LOGICAL                                            :: my_kp, my_rtp, my_xas
     436             :       TYPE(section_vals_type), POINTER                   :: print_key
     437             : 
     438         592 :       my_kp = .FALSE.
     439         592 :       my_xas = .FALSE.
     440         592 :       my_rtp = .FALSE.
     441         592 :       IF (PRESENT(kp)) my_kp = kp
     442         592 :       IF (PRESENT(xas)) my_xas = xas
     443         592 :       IF (PRESENT(rtp)) my_rtp = rtp
     444             : 
     445         592 :       exist = .FALSE.
     446         592 :       CALL section_vals_val_get(section, "WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
     447         592 :       IF (n_rep_val > 0) THEN
     448         462 :          CALL section_vals_val_get(section, "WFN_RESTART_FILE_NAME", c_val=filename)
     449             :       ELSE
     450         130 :          IF (my_xas) THEN
     451             :             ! try to read from the filename that is generated automatically from the printkey
     452           4 :             print_key => section_vals_get_subs_vals(section, "PRINT%RESTART")
     453             :             filename = cp_print_key_generate_filename(logger, print_key, &
     454           4 :                                                       extension="", my_local=.FALSE.)
     455         126 :          ELSE IF (my_rtp) THEN
     456             :             ! try to read from the filename that is generated automatically from the printkey
     457           3 :             print_key => section_vals_get_subs_vals(section, "REAL_TIME_PROPAGATION%PRINT%RESTART")
     458             :             filename = cp_print_key_generate_filename(logger, print_key, &
     459           3 :                                                       extension=".rtpwfn", my_local=.FALSE.)
     460         123 :          ELSE IF (my_kp) THEN
     461             :             ! try to read from the filename that is generated automatically from the printkey
     462           7 :             print_key => section_vals_get_subs_vals(section, "SCF%PRINT%RESTART")
     463             :             filename = cp_print_key_generate_filename(logger, print_key, &
     464           7 :                                                       extension=".kp", my_local=.FALSE.)
     465             :          ELSE
     466             :             ! try to read from the filename that is generated automatically from the printkey
     467         116 :             print_key => section_vals_get_subs_vals(section, "SCF%PRINT%RESTART")
     468             :             filename = cp_print_key_generate_filename(logger, print_key, &
     469         116 :                                                       extension=".wfn", my_local=.FALSE.)
     470             :          END IF
     471             :       END IF
     472         592 :       IF (.NOT. my_xas) THEN
     473         586 :          INQUIRE (FILE=filename, exist=exist)
     474             :       END IF
     475             : 
     476         592 :    END SUBROUTINE wfn_restart_file_name
     477             : 
     478             : ! **************************************************************************************************
     479             : !> \brief ...
     480             : !> \param mo_array ...
     481             : !> \param atomic_kind_set ...
     482             : !> \param qs_kind_set ...
     483             : !> \param particle_set ...
     484             : !> \param para_env ...
     485             : !> \param id_nr ...
     486             : !> \param multiplicity ...
     487             : !> \param dft_section ...
     488             : !> \param natom_mismatch ...
     489             : !> \param cdft ...
     490             : !> \param out_unit ...
     491             : ! **************************************************************************************************
     492         519 :    SUBROUTINE read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, &
     493             :                                        para_env, id_nr, multiplicity, dft_section, natom_mismatch, &
     494             :                                        cdft, out_unit)
     495             : 
     496             :       TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT)     :: mo_array
     497             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     498             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     499             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     500             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     501             :       INTEGER, INTENT(IN)                                :: id_nr, multiplicity
     502             :       TYPE(section_vals_type), POINTER                   :: dft_section
     503             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: natom_mismatch
     504             :       LOGICAL, INTENT(IN), OPTIONAL                      :: cdft
     505             :       INTEGER, INTENT(IN), OPTIONAL                      :: out_unit
     506             : 
     507             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_mo_set_from_restart'
     508             : 
     509             :       CHARACTER(LEN=default_path_length)                 :: file_name
     510             :       INTEGER                                            :: handle, ispin, my_out_unit, natom, &
     511             :                                                             nspin, restart_unit
     512             :       LOGICAL                                            :: exist, my_cdft
     513             :       TYPE(cp_logger_type), POINTER                      :: logger
     514             : 
     515         519 :       CALL timeset(routineN, handle)
     516         519 :       logger => cp_get_default_logger()
     517         519 :       my_cdft = .FALSE.
     518         519 :       IF (PRESENT(cdft)) my_cdft = cdft
     519         519 :       my_out_unit = -1
     520         519 :       IF (PRESENT(out_unit)) my_out_unit = out_unit
     521             : 
     522         519 :       nspin = SIZE(mo_array)
     523         519 :       restart_unit = -1
     524             : 
     525         519 :       IF (para_env%is_source()) THEN
     526             : 
     527         278 :          natom = SIZE(particle_set, 1)
     528         278 :          CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
     529         278 :          IF (id_nr /= 0) THEN
     530             :             ! Is it one of the backup files?
     531           1 :             file_name = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(id_nr))
     532             :          END IF
     533             : 
     534             :          CALL open_file(file_name=file_name, &
     535             :                         file_action="READ", &
     536             :                         file_form="UNFORMATTED", &
     537             :                         file_status="OLD", &
     538         278 :                         unit_number=restart_unit)
     539             : 
     540             :       END IF
     541             : 
     542             :       CALL read_mos_restart_low(mo_array, para_env=para_env, qs_kind_set=qs_kind_set, &
     543             :                                 particle_set=particle_set, natom=natom, &
     544         519 :                                 rst_unit=restart_unit, multiplicity=multiplicity, natom_mismatch=natom_mismatch)
     545             : 
     546         519 :       IF (PRESENT(natom_mismatch)) THEN
     547             :          ! read_mos_restart_low only the io_node returns natom_mismatch, must broadcast it
     548         489 :          CALL para_env%bcast(natom_mismatch)
     549         489 :          IF (natom_mismatch) THEN
     550           0 :             IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
     551           0 :             CALL timestop(handle)
     552           0 :             RETURN
     553             :          END IF
     554             :       END IF
     555             : 
     556             :       ! Close restart file
     557         519 :       IF (para_env%is_source()) THEN
     558         278 :          IF (my_out_unit > 0) THEN
     559             :             WRITE (UNIT=my_out_unit, FMT="(T2,A)") &
     560           6 :                "WFN_RESTART| Restart file "//TRIM(file_name)//" read"
     561             :          END IF
     562         278 :          CALL close_file(unit_number=restart_unit)
     563             :       END IF
     564             : 
     565             :       ! CDFT has no real dft_section and does not need to print
     566         519 :       IF (.NOT. my_cdft) THEN
     567        1394 :          DO ispin = 1, nspin
     568             :             CALL write_mo_set_to_output_unit(mo_array(ispin), atomic_kind_set, qs_kind_set, &
     569        1394 :                                              particle_set, dft_section, 4, 0, final_mos=.FALSE.)
     570             :          END DO
     571             :       END IF
     572             : 
     573         519 :       CALL timestop(handle)
     574             : 
     575             :    END SUBROUTINE read_mo_set_from_restart
     576             : 
     577             : ! **************************************************************************************************
     578             : !> \brief ...
     579             : !> \param mo_array ...
     580             : !> \param rt_mos ...
     581             : !> \param atomic_kind_set ...
     582             : !> \param qs_kind_set ...
     583             : !> \param particle_set ...
     584             : !> \param para_env ...
     585             : !> \param id_nr ...
     586             : !> \param multiplicity ...
     587             : !> \param dft_section ...
     588             : ! **************************************************************************************************
     589           8 :    SUBROUTINE read_rt_mos_from_restart(mo_array, rt_mos, atomic_kind_set, qs_kind_set, &
     590             :                                        particle_set, para_env, id_nr, multiplicity, dft_section)
     591             : 
     592             :       TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT)     :: mo_array
     593             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: rt_mos
     594             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     595             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     596             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     597             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     598             :       INTEGER, INTENT(IN)                                :: id_nr, multiplicity
     599             :       TYPE(section_vals_type), POINTER                   :: dft_section
     600             : 
     601             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_rt_mos_from_restart'
     602             : 
     603             :       CHARACTER(LEN=default_path_length)                 :: file_name
     604             :       INTEGER                                            :: handle, ispin, natom, nspin, &
     605             :                                                             restart_unit, unit_nr
     606             :       LOGICAL                                            :: exist
     607             :       TYPE(cp_logger_type), POINTER                      :: logger
     608             : 
     609           8 :       CALL timeset(routineN, handle)
     610           8 :       logger => cp_get_default_logger()
     611             : 
     612           8 :       nspin = SIZE(mo_array)
     613           8 :       restart_unit = -1
     614             : 
     615           8 :       IF (para_env%is_source()) THEN
     616             : 
     617           4 :          natom = SIZE(particle_set, 1)
     618           4 :          CALL wfn_restart_file_name(file_name, exist, dft_section, logger, rtp=.TRUE.)
     619           4 :          IF (id_nr /= 0) THEN
     620             :             ! Is it one of the backup files?
     621           0 :             file_name = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(id_nr))
     622             :          END IF
     623             : 
     624           4 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     625           4 :          IF (unit_nr > 0) THEN
     626           4 :             WRITE (unit_nr, '(T2,A)') "Read RTP restart from the file: "//TRIM(file_name)
     627             :          END IF
     628             : 
     629             :          CALL open_file(file_name=file_name, &
     630             :                         file_action="READ", &
     631             :                         file_form="UNFORMATTED", &
     632             :                         file_status="OLD", &
     633           4 :                         unit_number=restart_unit)
     634             : 
     635             :       END IF
     636             : 
     637             :       CALL read_mos_restart_low(mo_array, rt_mos=rt_mos, para_env=para_env, &
     638             :                                 particle_set=particle_set, qs_kind_set=qs_kind_set, natom=natom, &
     639           8 :                                 rst_unit=restart_unit, multiplicity=multiplicity)
     640             : 
     641             :       ! Close restart file
     642           8 :       IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
     643             : 
     644          16 :       DO ispin = 1, nspin
     645             :          CALL write_mo_set_to_output_unit(mo_array(ispin), atomic_kind_set, qs_kind_set, &
     646          16 :                                           particle_set, dft_section, 4, 0, final_mos=.FALSE.)
     647             :       END DO
     648             : 
     649           8 :       CALL timestop(handle)
     650             : 
     651           8 :    END SUBROUTINE read_rt_mos_from_restart
     652             : 
     653             : ! **************************************************************************************************
     654             : !> \brief Reading the mos from apreviously defined restart file
     655             : !> \param mos ...
     656             : !> \param para_env ...
     657             : !> \param qs_kind_set ...
     658             : !> \param particle_set ...
     659             : !> \param natom ...
     660             : !> \param rst_unit ...
     661             : !> \param multiplicity ...
     662             : !> \param rt_mos ...
     663             : !> \param natom_mismatch ...
     664             : !> \par History
     665             : !>      12.2007 created [MI]
     666             : !> \author MI
     667             : ! **************************************************************************************************
     668         583 :    SUBROUTINE read_mos_restart_low(mos, para_env, qs_kind_set, particle_set, natom, rst_unit, &
     669             :                                    multiplicity, rt_mos, natom_mismatch)
     670             : 
     671             :       TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT)     :: mos
     672             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     673             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     674             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     675             :       INTEGER, INTENT(IN)                                :: natom, rst_unit
     676             :       INTEGER, INTENT(in), OPTIONAL                      :: multiplicity
     677             :       TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER  :: rt_mos
     678             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: natom_mismatch
     679             : 
     680             :       INTEGER :: homo, homo_read, i, iatom, ikind, imat, irow, iset, iset_read, ishell, &
     681             :          ishell_read, iso, ispin, lfomo_read, lmax, lshell, my_mult, nao, nao_read, natom_read, &
     682             :          nelectron, nelectron_read, nmo, nmo_read, nnshell, nset, nset_max, nshell_max, nspin, &
     683             :          nspin_read, offset_read
     684         583 :       INTEGER, DIMENSION(:), POINTER                     :: nset_info, nshell
     685         583 :       INTEGER, DIMENSION(:, :), POINTER                  :: l, nshell_info
     686         583 :       INTEGER, DIMENSION(:, :, :), POINTER               :: nso_info, offset_info
     687             :       LOGICAL                                            :: minbas, natom_match, use_this
     688         583 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eig_read, occ_read
     689         583 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer, vecbuffer_read
     690             :       TYPE(cp_logger_type), POINTER                      :: logger
     691             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     692             :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
     693             : 
     694        1166 :       logger => cp_get_default_logger()
     695             : 
     696         583 :       nspin = SIZE(mos)
     697         583 :       nao = mos(1)%nao
     698         583 :       my_mult = 0
     699         583 :       IF (PRESENT(multiplicity)) my_mult = multiplicity
     700             : 
     701         583 :       IF (para_env%is_source()) THEN
     702         310 :          READ (rst_unit) natom_read, nspin_read, nao_read, nset_max, nshell_max
     703         310 :          IF (PRESENT(rt_mos)) THEN
     704           4 :             IF (nspin_read /= nspin) THEN
     705           0 :                CPABORT("To change nspin is not possible. ")
     706             :             END IF
     707             :          ELSE
     708             :             ! we should allow for restarting with different spin settings
     709         306 :             IF (nspin_read /= nspin) THEN
     710             :                WRITE (cp_logger_get_default_unit_nr(logger), *) &
     711           0 :                   "READ RESTART : WARNING : nspin is not equal "
     712             :             END IF
     713             :             ! this case needs fixing of homo/lfomo/nelec/occupations ...
     714         306 :             IF (nspin_read > nspin) THEN
     715           0 :                CPABORT("Reducing nspin is not possible. ")
     716             :             END IF
     717             :          END IF
     718             : 
     719         310 :          natom_match = (natom_read == natom)
     720             : 
     721         310 :          IF (natom_match) THEN ! actually do the read read
     722             : 
     723             :             ! Let's make it possible to change the basis set
     724        1550 :             ALLOCATE (nso_info(nshell_max, nset_max, natom_read))
     725        1240 :             ALLOCATE (nshell_info(nset_max, natom_read))
     726         930 :             ALLOCATE (nset_info(natom_read))
     727        1240 :             ALLOCATE (offset_info(nshell_max, nset_max, natom_read))
     728             : 
     729         310 :             IF (nao_read /= nao) THEN
     730             :                WRITE (cp_logger_get_default_unit_nr(logger), *) &
     731           1 :                   " READ RESTART : WARNING : DIFFERENT # AOs ", nao, nao_read
     732           1 :                IF (PRESENT(rt_mos)) &
     733           0 :                   CPABORT("To change basis is not possible. ")
     734             :             END IF
     735             : 
     736        1136 :             READ (rst_unit) nset_info
     737        2855 :             READ (rst_unit) nshell_info
     738        6900 :             READ (rst_unit) nso_info
     739             : 
     740         310 :             i = 1
     741        1136 :             DO iatom = 1, natom
     742        2671 :                DO iset = 1, nset_info(iatom)
     743        4878 :                   DO ishell = 1, nshell_info(iset, iatom)
     744        2517 :                      offset_info(ishell, iset, iatom) = i
     745        4052 :                      i = i + nso_info(ishell, iset, iatom)
     746             :                   END DO
     747             :                END DO
     748             :             END DO
     749             : 
     750         930 :             ALLOCATE (vecbuffer_read(1, nao_read))
     751             : 
     752             :          END IF ! natom_match
     753             :       END IF ! ionode
     754             : 
     755             :       ! make natom_match and natom_mismatch uniform across all nodes
     756         583 :       CALL para_env%bcast(natom_match)
     757         583 :       IF (PRESENT(natom_mismatch)) natom_mismatch = .NOT. natom_match
     758             :       ! handle natom_match false
     759         583 :       IF (.NOT. natom_match) THEN
     760           0 :          IF (PRESENT(natom_mismatch)) THEN
     761             :             WRITE (cp_logger_get_default_unit_nr(logger), *) &
     762           0 :                " READ RESTART : WARNING : DIFFERENT natom, returning ", natom, natom_read
     763             :             RETURN
     764             :          ELSE
     765           0 :             CPABORT("Incorrect number of atoms in restart file. ")
     766             :          END IF
     767             :       END IF
     768             : 
     769         583 :       CALL para_env%bcast(nspin_read)
     770             : 
     771        1749 :       ALLOCATE (vecbuffer(1, nao))
     772             : 
     773        1600 :       DO ispin = 1, nspin
     774             : 
     775        1017 :          nmo = mos(ispin)%nmo
     776        1017 :          homo = mos(ispin)%homo
     777        5330 :          mos(ispin)%eigenvalues(:) = 0.0_dp
     778        5330 :          mos(ispin)%occupation_numbers(:) = 0.0_dp
     779        1017 :          CALL cp_fm_set_all(mos(ispin)%mo_coeff, 0.0_dp)
     780             : 
     781        1017 :          IF (para_env%is_source() .AND. (nmo > 0)) THEN
     782         536 :             READ (rst_unit) nmo_read, homo_read, lfomo_read, nelectron_read
     783        2144 :             ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
     784        2744 :             eig_read = 0.0_dp
     785        2744 :             occ_read = 0.0_dp
     786             : 
     787         536 :             nmo = MIN(nmo, nmo_read)
     788             :             IF (nmo_read < nmo) &
     789             :                CALL cp_warn(__LOCATION__, &
     790             :                             "The number of MOs on the restart unit is smaller than the number of "// &
     791             :                             "the allocated MOs. The MO set will be padded with zeros!")
     792         536 :             IF (nmo_read > nmo) &
     793             :                CALL cp_warn(__LOCATION__, &
     794             :                             "The number of MOs on the restart unit is greater than the number of "// &
     795           7 :                             "the allocated MOs. The read MO set will be truncated!")
     796             : 
     797         536 :             READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read)
     798        2720 :             mos(ispin)%eigenvalues(1:nmo) = eig_read(1:nmo)
     799        2720 :             mos(ispin)%occupation_numbers(1:nmo) = occ_read(1:nmo)
     800         536 :             DEALLOCATE (eig_read, occ_read)
     801             : 
     802         536 :             mos(ispin)%homo = homo_read
     803         536 :             mos(ispin)%lfomo = lfomo_read
     804         536 :             IF (homo_read > nmo) THEN
     805           0 :                IF (nelectron_read == mos(ispin)%nelectron) THEN
     806             :                   CALL cp_warn(__LOCATION__, &
     807             :                                "The number of occupied MOs on the restart unit is larger than "// &
     808           0 :                                "the allocated MOs. The read MO set will be truncated and the occupation numbers recalculated!")
     809           0 :                   CALL set_mo_occupation(mo_set=mos(ispin))
     810             :                ELSE
     811             :                   ! can not make this a warning i.e. homo must be smaller than nmo
     812             :                   ! otherwise e.g. set_mo_occupation will go out of bounds
     813           0 :                   CPABORT("Number of occupied MOs on restart unit larger than allocated MOs. ")
     814             :                END IF
     815             :             END IF
     816             :          END IF
     817             : 
     818        1017 :          CALL para_env%bcast(nmo)
     819        1017 :          CALL para_env%bcast(mos(ispin)%homo)
     820        1017 :          CALL para_env%bcast(mos(ispin)%lfomo)
     821        1017 :          CALL para_env%bcast(mos(ispin)%nelectron)
     822        9643 :          CALL para_env%bcast(mos(ispin)%eigenvalues)
     823        9643 :          CALL para_env%bcast(mos(ispin)%occupation_numbers)
     824             : 
     825        1017 :          IF (PRESENT(rt_mos)) THEN
     826          24 :             DO imat = 2*ispin - 1, 2*ispin
     827          40 :                DO i = 1, nmo
     828          16 :                   IF (para_env%is_source()) THEN
     829         168 :                      READ (rst_unit) vecbuffer
     830             :                   ELSE
     831          88 :                      vecbuffer(1, :) = 0.0_dp
     832             :                   END IF
     833         656 :                   CALL para_env%bcast(vecbuffer)
     834             :                   CALL cp_fm_set_submatrix(rt_mos(imat), &
     835          32 :                                            vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
     836             :                END DO
     837             :             END DO
     838             :          ELSE
     839        5278 :             DO i = 1, nmo
     840        4269 :                IF (para_env%is_source()) THEN
     841      145910 :                   READ (rst_unit) vecbuffer_read
     842             :                   ! now, try to assign the read to the real vector
     843             :                   ! in case the basis set changed this involves some guessing
     844        2180 :                   irow = 1
     845        9754 :                   DO iatom = 1, natom
     846        7574 :                      NULLIFY (orb_basis_set, dftb_parameter, l, nshell)
     847        7574 :                      CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     848             :                      CALL get_qs_kind(qs_kind_set(ikind), &
     849        7574 :                                       basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
     850        7574 :                      IF (ASSOCIATED(orb_basis_set)) THEN
     851             :                         CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     852             :                                                nset=nset, &
     853             :                                                nshell=nshell, &
     854        7526 :                                                l=l)
     855        7526 :                         minbas = .FALSE.
     856          48 :                      ELSEIF (ASSOCIATED(dftb_parameter)) THEN
     857          48 :                         CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
     858          48 :                         nset = 1
     859          48 :                         minbas = .TRUE.
     860             :                      ELSE
     861             :                         ! assume an atom without basis set
     862             :                         ! CPABORT("Unknown basis set type. ")
     863           0 :                         nset = 0
     864             :                      END IF
     865             : 
     866        7574 :                      use_this = .TRUE.
     867        7574 :                      iset_read = 1
     868       35460 :                      DO iset = 1, nset
     869       18132 :                         ishell_read = 1
     870       18132 :                         IF (minbas) THEN
     871          48 :                            nnshell = lmax + 1
     872             :                         ELSE
     873       18084 :                            nnshell = nshell(iset)
     874             :                         END IF
     875       57349 :                         DO ishell = 1, nnshell
     876       31643 :                            IF (minbas) THEN
     877          72 :                               lshell = ishell - 1
     878             :                            ELSE
     879       31571 :                               lshell = l(ishell, iset)
     880             :                            END IF
     881       31643 :                            IF (iset_read > nset_info(iatom)) use_this = .FALSE.
     882             :                            IF (use_this) THEN ! avoids out of bound access of the lower line if false
     883       31619 :                               IF (nso(lshell) == nso_info(ishell_read, iset_read, iatom)) THEN
     884       31619 :                                  offset_read = offset_info(ishell_read, iset_read, iatom)
     885       31619 :                                  ishell_read = ishell_read + 1
     886       31619 :                                  IF (ishell_read > nshell_info(iset, iatom)) THEN
     887       18128 :                                     ishell_read = 1
     888       18128 :                                     iset_read = iset_read + 1
     889             :                                  END IF
     890             :                               ELSE
     891             :                                  use_this = .FALSE.
     892             :                               END IF
     893             :                            END IF
     894      103652 :                            DO iso = 1, nso(lshell)
     895       72009 :                               IF (use_this) THEN
     896       71865 :                                  IF (offset_read - 1 + iso .LT. 1 .OR. offset_read - 1 + iso .GT. nao_read) THEN
     897           0 :                                     vecbuffer(1, irow) = 0.0_dp
     898             :                                  ELSE
     899       71865 :                                     vecbuffer(1, irow) = vecbuffer_read(1, offset_read - 1 + iso)
     900             :                                  END IF
     901             :                               ELSE
     902         144 :                                  vecbuffer(1, irow) = 0.0_dp
     903             :                               END IF
     904      103652 :                               irow = irow + 1
     905             :                            END DO
     906       49775 :                            use_this = .TRUE.
     907             :                         END DO
     908             :                      END DO
     909             :                   END DO
     910             : 
     911             :                ELSE
     912             : 
     913       72620 :                   vecbuffer(1, :) = 0.0_dp
     914             : 
     915             :                END IF
     916             : 
     917      574429 :                CALL para_env%bcast(vecbuffer)
     918             :                CALL cp_fm_set_submatrix(mos(ispin)%mo_coeff, &
     919        5278 :                                         vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
     920             :             END DO
     921             :          END IF
     922             :          ! Skip extra MOs if there any
     923        1017 :          IF (para_env%is_source()) THEN
     924             :             !ignore nmo = 0
     925         539 :             IF (nmo > 0) THEN
     926         560 :                DO i = nmo + 1, nmo_read
     927        1504 :                   READ (rst_unit) vecbuffer_read
     928             :                END DO
     929             :             END IF
     930             :          END IF
     931             : 
     932        1600 :          IF (.NOT. PRESENT(rt_mos)) THEN
     933        1009 :             IF (ispin == 1 .AND. nspin_read < nspin) THEN
     934             : 
     935           0 :                mos(ispin + 1)%homo = mos(ispin)%homo
     936           0 :                mos(ispin + 1)%lfomo = mos(ispin)%lfomo
     937           0 :                nelectron = mos(ispin)%nelectron
     938           0 :                IF (my_mult .NE. 1) THEN
     939             :                   CALL cp_abort(__LOCATION__, &
     940           0 :                                 "Restarting an LSD calculation from an LDA wfn only works for multiplicity=1 (singlets).")
     941             :                END IF
     942           0 :                IF (mos(ispin + 1)%nelectron < 0) THEN
     943           0 :                   CPABORT("LSD: too few electrons for this multiplisity. ")
     944             :                END IF
     945           0 :                mos(ispin + 1)%eigenvalues = mos(ispin)%eigenvalues
     946           0 :                mos(ispin)%occupation_numbers = mos(ispin)%occupation_numbers/2.0_dp
     947           0 :                mos(ispin + 1)%occupation_numbers = mos(ispin)%occupation_numbers
     948           0 :                CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos(ispin + 1)%mo_coeff)
     949           0 :                EXIT
     950             :             END IF
     951             :          END IF
     952             :       END DO ! ispin
     953             : 
     954         583 :       DEALLOCATE (vecbuffer)
     955             : 
     956         583 :       IF (para_env%is_source()) THEN
     957         310 :          DEALLOCATE (vecbuffer_read)
     958         310 :          DEALLOCATE (offset_info)
     959         310 :          DEALLOCATE (nso_info)
     960         310 :          DEALLOCATE (nshell_info)
     961         310 :          DEALLOCATE (nset_info)
     962             :       END IF
     963             : 
     964        1166 :    END SUBROUTINE read_mos_restart_low
     965             : 
     966             : ! **************************************************************************************************
     967             : !> \brief Write MO information to output file (eigenvalues, occupation numbers, coefficients)
     968             : !> \param mo_set ...
     969             : !> \param atomic_kind_set ...
     970             : !> \param qs_kind_set ...
     971             : !> \param particle_set ...
     972             : !> \param dft_section ...
     973             : !> \param before Digits before the dot
     974             : !> \param kpoint An integer that labels the current k point, e.g. its index
     975             : !> \param final_mos ...
     976             : !> \param spin ...
     977             : !> \param solver_method ...
     978             : !> \param rtp ...
     979             : !> \param cpart ...
     980             : !> \param sim_step ...
     981             : !> \param umo_set ...
     982             : !> \date    15.05.2001
     983             : !> \par History:
     984             : !>       - Optionally print Cartesian MOs (20.04.2005, MK)
     985             : !>       - Revise printout of MO information (05.05.2021, MK)
     986             : !> \par Variables
     987             : !>       - after : Number of digits after point.
     988             : !>       - before: Number of digits before point.
     989             : !> \author  Matthias Krack (MK)
     990             : !> \version 1.1
     991             : ! **************************************************************************************************
     992        7679 :    SUBROUTINE write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
     993             :                                           dft_section, before, kpoint, final_mos, spin, &
     994             :                                           solver_method, rtp, cpart, sim_step, umo_set)
     995             : 
     996             :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     997             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     998             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     999             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1000             :       TYPE(section_vals_type), POINTER                   :: dft_section
    1001             :       INTEGER, INTENT(IN)                                :: before, kpoint
    1002             :       LOGICAL, INTENT(IN), OPTIONAL                      :: final_mos
    1003             :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: spin
    1004             :       CHARACTER(LEN=2), INTENT(IN), OPTIONAL             :: solver_method
    1005             :       LOGICAL, INTENT(IN), OPTIONAL                      :: rtp
    1006             :       INTEGER, INTENT(IN), OPTIONAL                      :: cpart, sim_step
    1007             :       TYPE(mo_set_type), INTENT(IN), OPTIONAL            :: umo_set
    1008             : 
    1009             :       CHARACTER(LEN=12)                                  :: symbol
    1010        7679 :       CHARACTER(LEN=12), DIMENSION(:), POINTER           :: bcgf_symbol
    1011             :       CHARACTER(LEN=14)                                  :: fmtstr5
    1012             :       CHARACTER(LEN=15)                                  :: energy_str, orbital_str, step_string
    1013             :       CHARACTER(LEN=2)                                   :: element_symbol, my_solver_method
    1014             :       CHARACTER(LEN=2*default_string_length)             :: name
    1015             :       CHARACTER(LEN=21)                                  :: vector_str
    1016             :       CHARACTER(LEN=22)                                  :: fmtstr4
    1017             :       CHARACTER(LEN=24)                                  :: fmtstr2
    1018             :       CHARACTER(LEN=25)                                  :: fmtstr1
    1019             :       CHARACTER(LEN=29)                                  :: fmtstr6
    1020             :       CHARACTER(LEN=4)                                   :: reim
    1021             :       CHARACTER(LEN=40)                                  :: fmtstr3
    1022        7679 :       CHARACTER(LEN=6), DIMENSION(:), POINTER            :: bsgf_symbol
    1023             :       INTEGER :: after, first_mo, from, homo, iatom, icgf, ico, icol, ikind, imo, irow, iset, &
    1024             :          isgf, ishell, iso, iw, jcol, last_mo, left, lmax, lshell, nao, natom, ncgf, ncol, nmo, &
    1025             :          nset, nsgf, numo, right, scf_step, to, width
    1026        7679 :       INTEGER, DIMENSION(:), POINTER                     :: mo_index_range, nshell
    1027        7679 :       INTEGER, DIMENSION(:, :), POINTER                  :: l
    1028             :       LOGICAL                                            :: ionode, my_final, my_rtp, &
    1029             :                                                             print_cartesian, print_eigvals, &
    1030             :                                                             print_eigvecs, print_occup, &
    1031             :                                                             should_output
    1032             :       REAL(KIND=dp)                                      :: gap, maxocc
    1033        7679 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: mo_eigenvalues, mo_occupation_numbers
    1034        7679 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cmatrix, smatrix
    1035        7679 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eigenvalues, occupation_numbers
    1036             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff, umo_coeff
    1037             :       TYPE(cp_logger_type), POINTER                      :: logger
    1038             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1039             :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
    1040             : 
    1041        7679 :       NULLIFY (bcgf_symbol)
    1042        7679 :       NULLIFY (bsgf_symbol)
    1043        7679 :       NULLIFY (logger)
    1044        7679 :       NULLIFY (mo_index_range)
    1045        7679 :       NULLIFY (nshell)
    1046        7679 :       NULLIFY (mo_coeff)
    1047             : 
    1048       15358 :       logger => cp_get_default_logger()
    1049        7679 :       ionode = logger%para_env%is_source()
    1050        7679 :       CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVALUES", l_val=print_eigvals)
    1051        7679 :       CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVECTORS", l_val=print_eigvecs)
    1052        7679 :       CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occup)
    1053        7679 :       CALL section_vals_val_get(dft_section, "PRINT%MO%CARTESIAN", l_val=print_cartesian)
    1054        7679 :       CALL section_vals_val_get(dft_section, "PRINT%MO%MO_INDEX_RANGE", i_vals=mo_index_range)
    1055        7679 :       CALL section_vals_val_get(dft_section, "PRINT%MO%NDIGITS", i_val=after)
    1056        7679 :       after = MIN(MAX(after, 1), 16)
    1057             : 
    1058             :       ! Do we print the final MO information after SCF convergence is reached (default: no)
    1059        7679 :       IF (PRESENT(final_mos)) THEN
    1060        7671 :          my_final = final_mos
    1061             :       ELSE
    1062             :          my_final = .FALSE.
    1063             :       END IF
    1064             : 
    1065             :       ! complex MOS for RTP, no eigenvalues
    1066        7679 :       my_rtp = .FALSE.
    1067        7679 :       IF (PRESENT(rtp)) THEN
    1068           8 :          my_rtp = rtp
    1069             :          ! print the first time step if MO print required
    1070             :          should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
    1071             :                                                           "PRINT%MO"), cp_p_file) &
    1072           8 :                          .OR. (sim_step == 1)
    1073             :       ELSE
    1074             :          should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
    1075        8592 :                                                           "PRINT%MO"), cp_p_file) .OR. my_final
    1076             :       END IF
    1077             : 
    1078        7679 :       IF ((.NOT. should_output) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup))) RETURN
    1079             : 
    1080        6756 :       IF (my_rtp) THEN
    1081           8 :          CPASSERT(PRESENT(sim_step))
    1082           8 :          CPASSERT(PRESENT(cpart))
    1083           8 :          scf_step = sim_step
    1084           8 :          IF (cpart == 0) THEN
    1085           4 :             reim = "IMAG"
    1086             :          ELSE
    1087           4 :             reim = "REAL"
    1088             :          END IF
    1089           8 :          print_eigvals = .FALSE.
    1090             :       ELSE
    1091        6748 :          scf_step = MAX(0, logger%iter_info%iteration(logger%iter_info%n_rlevel) - 1)
    1092             :       END IF
    1093             : 
    1094        6756 :       IF (.NOT. my_final) THEN
    1095        5618 :          IF (.NOT. my_rtp) THEN
    1096        5610 :             step_string = " AFTER SCF STEP"
    1097             :          ELSE
    1098           8 :             step_string = " AFTER RTP STEP"
    1099             :          END IF
    1100             :       END IF
    1101             : 
    1102        6756 :       IF (PRESENT(solver_method)) THEN
    1103        6602 :          my_solver_method = solver_method
    1104             :       ELSE
    1105             :          ! Traditional diagonalization is assumed as default solver method
    1106         154 :          my_solver_method = "TD"
    1107             :       END IF
    1108             : 
    1109             :       ! Retrieve MO information
    1110             :       CALL get_mo_set(mo_set=mo_set, &
    1111             :                       mo_coeff=mo_coeff, &
    1112             :                       eigenvalues=eigenvalues, &
    1113             :                       occupation_numbers=occupation_numbers, &
    1114             :                       homo=homo, &
    1115             :                       maxocc=maxocc, &
    1116             :                       nao=nao, &
    1117        6756 :                       nmo=nmo)
    1118        6756 :       IF (PRESENT(umo_set)) THEN
    1119             :          CALL get_mo_set(mo_set=umo_set, &
    1120             :                          mo_coeff=umo_coeff, &
    1121          20 :                          nmo=numo)
    1122          20 :          nmo = nmo + numo
    1123             :       ELSE
    1124        6736 :          numo = 0
    1125             :       END IF
    1126       20218 :       ALLOCATE (mo_eigenvalues(nmo))
    1127       54320 :       mo_eigenvalues(:) = 0.0_dp
    1128       54210 :       mo_eigenvalues(1:nmo - numo) = eigenvalues(1:nmo - numo)
    1129       13462 :       ALLOCATE (mo_occupation_numbers(nmo))
    1130       54320 :       mo_occupation_numbers(:) = 0.0_dp
    1131       54210 :       mo_occupation_numbers(1:nmo - numo) = occupation_numbers(1:nmo - numo)
    1132        6756 :       IF (numo > 0) THEN
    1133             :          CALL get_mo_set(mo_set=umo_set, &
    1134          20 :                          eigenvalues=eigenvalues)
    1135         130 :          mo_eigenvalues(nmo - numo + 1:nmo) = eigenvalues(1:numo)
    1136             :       END IF
    1137             : 
    1138        6756 :       IF (print_eigvecs) THEN
    1139       20494 :          ALLOCATE (smatrix(nao, nmo))
    1140        5136 :          CALL cp_fm_get_submatrix(mo_coeff, smatrix(1:nao, 1:nmo - numo))
    1141        5136 :          IF (numo > 0) THEN
    1142          14 :             CALL cp_fm_get_submatrix(umo_coeff, smatrix(1:nao, nmo - numo + 1:nmo))
    1143             :          END IF
    1144        5136 :          IF (.NOT. ionode) THEN
    1145        2568 :             DEALLOCATE (smatrix)
    1146             :          END IF
    1147             :       END IF
    1148             : 
    1149             :       iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO", &
    1150             :                                 ignore_should_output=should_output, &
    1151        6756 :                                 extension=".MOLog")
    1152             : 
    1153        6756 :       IF (iw > 0) THEN
    1154             : 
    1155        3378 :          CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
    1156        3378 :          CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
    1157             : 
    1158             :          ! Definition of the variable formats
    1159             : 
    1160        3378 :          fmtstr1 = "(T2,A,21X,  (  X,I5,  X))"
    1161        3378 :          fmtstr2 = "(T2,A,21X,  (1X,F  .  ))"
    1162        3378 :          fmtstr3 = "(T2,A,I5,1X,I5,1X,A,1X,A6,  (1X,F  .  ))"
    1163             : 
    1164        3378 :          width = before + after + 3
    1165        3378 :          ncol = INT(56/width)
    1166             : 
    1167        3378 :          right = MAX((after - 2), 1)
    1168        3378 :          left = width - right - 5
    1169             : 
    1170        3378 :          WRITE (UNIT=fmtstr1(11:12), FMT="(I2)") ncol
    1171        3378 :          WRITE (UNIT=fmtstr1(14:15), FMT="(I2)") left
    1172        3378 :          WRITE (UNIT=fmtstr1(21:22), FMT="(I2)") right
    1173             : 
    1174        3378 :          WRITE (UNIT=fmtstr2(11:12), FMT="(I2)") ncol
    1175        3378 :          WRITE (UNIT=fmtstr2(18:19), FMT="(I2)") width - 1
    1176        3378 :          WRITE (UNIT=fmtstr2(21:22), FMT="(I2)") after
    1177             : 
    1178        3378 :          WRITE (UNIT=fmtstr3(27:28), FMT="(I2)") ncol
    1179        3378 :          WRITE (UNIT=fmtstr3(34:35), FMT="(I2)") width - 1
    1180        3378 :          WRITE (UNIT=fmtstr3(37:38), FMT="(I2)") after
    1181             : 
    1182        3378 :          IF (my_final .OR. (my_solver_method == "TD")) THEN
    1183        2551 :             energy_str = "EIGENVALUES"
    1184        2551 :             vector_str = "EIGENVECTORS"
    1185             :          ELSE
    1186         827 :             energy_str = "ENERGIES"
    1187         827 :             vector_str = "COEFFICIENTS"
    1188             :          END IF
    1189             : 
    1190        3378 :          IF (my_rtp) THEN
    1191           4 :             energy_str = "ZEROS"
    1192           4 :             vector_str = TRIM(reim)//" RTP COEFFICIENTS"
    1193             :          END IF
    1194             : 
    1195        3378 :          IF (print_eigvecs) THEN
    1196             : 
    1197        2568 :             IF (print_cartesian) THEN
    1198             : 
    1199          85 :                orbital_str = "CARTESIAN"
    1200             : 
    1201         340 :                ALLOCATE (cmatrix(ncgf, ncgf))
    1202       28495 :                cmatrix = 0.0_dp
    1203             : 
    1204             :                ! Transform spherical MOs to Cartesian MOs
    1205          85 :                icgf = 1
    1206          85 :                isgf = 1
    1207         271 :                DO iatom = 1, natom
    1208         186 :                   NULLIFY (orb_basis_set, dftb_parameter)
    1209         186 :                   CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
    1210             :                   CALL get_qs_kind(qs_kind_set(ikind), &
    1211             :                                    basis_set=orb_basis_set, &
    1212         186 :                                    dftb_parameter=dftb_parameter)
    1213         457 :                   IF (ASSOCIATED(orb_basis_set)) THEN
    1214             :                      CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1215             :                                             nset=nset, &
    1216             :                                             nshell=nshell, &
    1217         150 :                                             l=l)
    1218         466 :                      DO iset = 1, nset
    1219         944 :                         DO ishell = 1, nshell(iset)
    1220         478 :                            lshell = l(ishell, iset)
    1221             :                            CALL dgemm("T", "N", nco(lshell), nmo, nso(lshell), 1.0_dp, &
    1222             :                                       orbtramat(lshell)%c2s, nso(lshell), &
    1223             :                                       smatrix(isgf, 1), nsgf, 0.0_dp, &
    1224         478 :                                       cmatrix(icgf, 1), ncgf)
    1225         478 :                            icgf = icgf + nco(lshell)
    1226         794 :                            isgf = isgf + nso(lshell)
    1227             :                         END DO
    1228             :                      END DO
    1229          36 :                   ELSE IF (ASSOCIATED(dftb_parameter)) THEN
    1230          36 :                      CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
    1231          90 :                      DO ishell = 1, lmax + 1
    1232          54 :                         lshell = ishell - 1
    1233             :                         CALL dgemm("T", "N", nco(lshell), nsgf, nso(lshell), 1.0_dp, &
    1234             :                                    orbtramat(lshell)%c2s, nso(lshell), &
    1235             :                                    smatrix(isgf, 1), nsgf, 0.0_dp, &
    1236          54 :                                    cmatrix(icgf, 1), ncgf)
    1237          54 :                         icgf = icgf + nco(lshell)
    1238          90 :                         isgf = isgf + nso(lshell)
    1239             :                      END DO
    1240             :                   ELSE
    1241             :                      ! assume atom without basis set
    1242             :                      ! CPABORT("Unknown basis set type")
    1243             :                   END IF
    1244             :                END DO ! iatom
    1245             : 
    1246             :             ELSE
    1247             : 
    1248        2483 :                orbital_str = "SPHERICAL"
    1249             : 
    1250             :             END IF ! print_cartesian
    1251             : 
    1252             :             name = TRIM(energy_str)//", OCCUPATION NUMBERS, AND "// &
    1253        2568 :                    TRIM(orbital_str)//" "//TRIM(vector_str)
    1254             : 
    1255        2568 :             IF (.NOT. my_final) &
    1256        2224 :                WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
    1257             : 
    1258         810 :          ELSE IF (print_occup .OR. print_eigvals) THEN
    1259         810 :             name = TRIM(energy_str)//" AND OCCUPATION NUMBERS"
    1260             : 
    1261         810 :             IF (.NOT. my_final) &
    1262         585 :                WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
    1263             :          END IF ! print_eigvecs
    1264             : 
    1265             :          ! Print headline
    1266        3378 :          IF (PRESENT(spin) .AND. (kpoint > 0)) THEN
    1267             :             WRITE (UNIT=iw, FMT="(/,T2,A,I0)") &
    1268           0 :                "MO| "//TRIM(spin)//" "//TRIM(name)//" FOR K POINT ", kpoint
    1269             :          ELSE IF (PRESENT(spin)) THEN
    1270             :             WRITE (UNIT=iw, FMT="(/,T2,A)") &
    1271         324 :                "MO| "//TRIM(spin)//" "//TRIM(name)
    1272        3054 :          ELSE IF (kpoint > 0) THEN
    1273             :             WRITE (UNIT=iw, FMT="(/,T2,A,I0)") &
    1274         101 :                "MO| "//TRIM(name)//" FOR K POINT ", kpoint
    1275             :          ELSE
    1276             :             WRITE (UNIT=iw, FMT="(/,T2,A)") &
    1277        2953 :                "MO| "//TRIM(name)
    1278             :          END IF
    1279             : 
    1280             :          ! Check if only a subset of the MOs has to be printed
    1281             :          IF ((mo_index_range(1) > 0) .AND. &
    1282        3378 :              (mo_index_range(2) > 0) .AND. &
    1283             :              (mo_index_range(2) >= mo_index_range(1))) THEN
    1284         171 :             first_mo = MAX(1, mo_index_range(1))
    1285         171 :             last_mo = MIN(nmo, mo_index_range(2))
    1286             :          ELSE
    1287        3207 :             first_mo = 1
    1288        3207 :             last_mo = nmo
    1289             :          END IF
    1290             : 
    1291        3378 :          IF (print_eigvecs) THEN
    1292             : 
    1293             :             ! Print full MO information
    1294             : 
    1295        5582 :             DO icol = first_mo, last_mo, ncol
    1296             : 
    1297        3014 :                from = icol
    1298        3014 :                to = MIN((from + ncol - 1), last_mo)
    1299             : 
    1300        3014 :                WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1301             :                WRITE (UNIT=iw, FMT=fmtstr1) &
    1302       14230 :                   "MO|", (jcol, jcol=from, to)
    1303             :                WRITE (UNIT=iw, FMT=fmtstr2) &
    1304        3014 :                   "MO|", (mo_eigenvalues(jcol), jcol=from, to)
    1305        3014 :                WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1306             :                WRITE (UNIT=iw, FMT=fmtstr2) &
    1307        3014 :                   "MO|", (mo_occupation_numbers(jcol), jcol=from, to)
    1308        3014 :                WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1309             : 
    1310        3014 :                irow = 1
    1311             : 
    1312       11959 :                DO iatom = 1, natom
    1313             : 
    1314        6377 :                   IF (iatom /= 1) WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1315             : 
    1316        6377 :                   NULLIFY (orb_basis_set, dftb_parameter)
    1317             :                   CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
    1318        6377 :                                        element_symbol=element_symbol, kind_number=ikind)
    1319             :                   CALL get_qs_kind(qs_kind_set(ikind), &
    1320             :                                    basis_set=orb_basis_set, &
    1321        6377 :                                    dftb_parameter=dftb_parameter)
    1322             : 
    1323       15768 :                   IF (print_cartesian) THEN
    1324             : 
    1325         852 :                      IF (ASSOCIATED(orb_basis_set)) THEN
    1326             :                         CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1327             :                                                nset=nset, &
    1328             :                                                nshell=nshell, &
    1329             :                                                l=l, &
    1330         744 :                                                cgf_symbol=bcgf_symbol)
    1331             : 
    1332         744 :                         icgf = 1
    1333        2520 :                         DO iset = 1, nset
    1334        4944 :                            DO ishell = 1, nshell(iset)
    1335        2424 :                               lshell = l(ishell, iset)
    1336       10824 :                               DO ico = 1, nco(lshell)
    1337             :                                  WRITE (UNIT=iw, FMT=fmtstr3) &
    1338        6624 :                                     "MO|", irow, iatom, ADJUSTR(element_symbol), bcgf_symbol(icgf), &
    1339       13248 :                                     (cmatrix(irow, jcol), jcol=from, to)
    1340        6624 :                                  icgf = icgf + 1
    1341        9048 :                                  irow = irow + 1
    1342             :                               END DO
    1343             :                            END DO
    1344             :                         END DO
    1345         108 :                      ELSE IF (ASSOCIATED(dftb_parameter)) THEN
    1346         108 :                         CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
    1347         108 :                         icgf = 1
    1348         270 :                         DO ishell = 1, lmax + 1
    1349         162 :                            lshell = ishell - 1
    1350         540 :                            DO ico = 1, nco(lshell)
    1351         270 :                               symbol = cgf_symbol(1, indco(1:3, icgf))
    1352         270 :                               symbol(1:2) = "  "
    1353             :                               WRITE (UNIT=iw, FMT=fmtstr3) &
    1354         270 :                                  "MO|", irow, iatom, ADJUSTR(element_symbol), symbol, &
    1355         540 :                                  (cmatrix(irow, jcol), jcol=from, to)
    1356         270 :                               icgf = icgf + 1
    1357         432 :                               irow = irow + 1
    1358             :                            END DO
    1359             :                         END DO
    1360             :                      ELSE
    1361             :                         ! assume atom without basis set
    1362             :                         ! CPABORT("Unknown basis set type")
    1363             :                      END IF
    1364             : 
    1365             :                   ELSE
    1366             : 
    1367        5525 :                      IF (ASSOCIATED(orb_basis_set)) THEN
    1368             :                         CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
    1369             :                                                nset=nset, &
    1370             :                                                nshell=nshell, &
    1371             :                                                l=l, &
    1372        5525 :                                                sgf_symbol=bsgf_symbol)
    1373        5525 :                         isgf = 1
    1374       16025 :                         DO iset = 1, nset
    1375       28923 :                            DO ishell = 1, nshell(iset)
    1376       12898 :                               lshell = l(ishell, iset)
    1377       54224 :                               DO iso = 1, nso(lshell)
    1378             :                                  WRITE (UNIT=iw, FMT=fmtstr3) &
    1379       30826 :                                     "MO|", irow, iatom, ADJUSTR(element_symbol), bsgf_symbol(isgf), &
    1380       61652 :                                     (smatrix(irow, jcol), jcol=from, to)
    1381       30826 :                                  isgf = isgf + 1
    1382       43724 :                                  irow = irow + 1
    1383             :                               END DO
    1384             :                            END DO
    1385             :                         END DO
    1386           0 :                      ELSE IF (ASSOCIATED(dftb_parameter)) THEN
    1387           0 :                         CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
    1388           0 :                         isgf = 1
    1389           0 :                         DO ishell = 1, lmax + 1
    1390           0 :                            lshell = ishell - 1
    1391           0 :                            DO iso = 1, nso(lshell)
    1392           0 :                               symbol = sgf_symbol(1, lshell, -lshell + iso - 1)
    1393           0 :                               symbol(1:2) = "  "
    1394             :                               WRITE (UNIT=iw, FMT=fmtstr3) &
    1395           0 :                                  "MO|", irow, iatom, ADJUSTR(element_symbol), symbol, &
    1396           0 :                                  (smatrix(irow, jcol), jcol=from, to)
    1397           0 :                               isgf = isgf + 1
    1398           0 :                               irow = irow + 1
    1399             :                            END DO
    1400             :                         END DO
    1401             :                      ELSE
    1402             :                         ! assume atom without basis set
    1403             :                         ! CPABORT("Unknown basis set type")
    1404             :                      END IF
    1405             : 
    1406             :                   END IF ! print_cartesian
    1407             : 
    1408             :                END DO ! iatom
    1409             : 
    1410             :             END DO ! icol
    1411             : 
    1412        2568 :             WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1413             : 
    1414             :             ! Release work storage
    1415             : 
    1416        2568 :             DEALLOCATE (smatrix)
    1417        2568 :             IF (print_cartesian) THEN
    1418          85 :                DEALLOCATE (cmatrix)
    1419             :             END IF
    1420             : 
    1421         810 :          ELSE IF (print_occup .OR. print_eigvals) THEN
    1422             : 
    1423         810 :             WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1424         810 :             fmtstr4 = "(T2,A,I7,3(1X,F22.  ))"
    1425         810 :             WRITE (UNIT=fmtstr4(19:20), FMT="(I2)") after
    1426         810 :             IF (my_final .OR. (my_solver_method == "TD")) THEN
    1427             :                WRITE (UNIT=iw, FMT="(A)") &
    1428         810 :                   " MO|  Index      Eigenvalue [a.u.]        Eigenvalue [eV]             Occupation"
    1429             :             ELSE
    1430             :                WRITE (UNIT=iw, FMT="(A)") &
    1431           0 :                   " MO|  Index          Energy [a.u.]            Energy [eV]             Occupation"
    1432             :             END IF
    1433       13332 :             DO imo = first_mo, last_mo
    1434             :                WRITE (UNIT=iw, FMT=fmtstr4) &
    1435       12522 :                   "MO|", imo, mo_eigenvalues(imo), &
    1436       12522 :                   mo_eigenvalues(imo)*evolt, &
    1437       25854 :                   mo_occupation_numbers(imo)
    1438             :             END DO
    1439         810 :             fmtstr5 = "(A,T59,F22.  )"
    1440         810 :             WRITE (UNIT=fmtstr5(12:13), FMT="(I2)") after
    1441             :             WRITE (UNIT=iw, FMT=fmtstr5) &
    1442         810 :                " MO| Sum:", accurate_sum(mo_occupation_numbers(:))
    1443             : 
    1444             :          END IF ! print_eigvecs
    1445             : 
    1446        3378 :          IF (.NOT. my_rtp) THEN
    1447        3374 :             fmtstr6 = "(A,T18,F17.  ,A,T41,F17.  ,A)"
    1448        3374 :             WRITE (UNIT=fmtstr6(12:13), FMT="(I2)") after
    1449        3374 :             WRITE (UNIT=fmtstr6(25:26), FMT="(I2)") after
    1450             :             WRITE (UNIT=iw, FMT=fmtstr6) &
    1451        3374 :                " MO| E(Fermi):", mo_set%mu, " a.u.", mo_set%mu*evolt, " eV"
    1452             :          END IF
    1453        3378 :          IF ((homo > 0) .AND. .NOT. my_rtp) THEN
    1454        3349 :             IF ((mo_occupation_numbers(homo) == maxocc) .AND. (last_mo > homo)) THEN
    1455             :                gap = mo_eigenvalues(homo + 1) - &
    1456         384 :                      mo_eigenvalues(homo)
    1457             :                WRITE (UNIT=iw, FMT=fmtstr6) &
    1458         384 :                   " MO| Band gap:", gap, " a.u.", gap*evolt, " eV"
    1459             :             END IF
    1460             :          END IF
    1461        3378 :          WRITE (UNIT=iw, FMT="(A)") ""
    1462             : 
    1463             :       END IF ! iw
    1464             : 
    1465        6756 :       IF (ALLOCATED(mo_eigenvalues)) DEALLOCATE (mo_eigenvalues)
    1466        6756 :       IF (ALLOCATED(mo_occupation_numbers)) DEALLOCATE (mo_occupation_numbers)
    1467             : 
    1468             :       CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%MO", &
    1469        6756 :                                         ignore_should_output=should_output)
    1470             : 
    1471       29793 :    END SUBROUTINE write_mo_set_to_output_unit
    1472             : 
    1473             : END MODULE qs_mo_io

Generated by: LCOV version 1.15