LCOV - code coverage report
Current view: top level - src - qs_mo_io.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 90.4 % 594 537
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 9 9

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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_create,&
      30              :                                               dbcsr_p_type,&
      31              :                                               dbcsr_release,&
      32              :                                               dbcsr_type
      33              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_checksum
      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       161937 :    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       161937 :       CALL timeset(routineN, handle)
     121       161937 :       logger => cp_get_default_logger()
     122              : 
     123              :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     124       161937 :                                            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        17325 :          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        12072 :             DO ispin = 1, SIZE(mo_array)
     132         6543 :                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        12072 :                                      mo_array(ispin)%mo_coeff) !fm->dbcsr
     137              :             END DO
     138              :          END IF
     139              : 
     140        51975 :          DO ikey = 1, SIZE(keys)
     141              : 
     142        34650 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     143       161937 :                                                  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        17337 :                                            do_backup=.TRUE., file_form="UNFORMATTED")
     148              : 
     149              :                CALL write_mo_set_low(mo_array, particle_set=particle_set, &
     150        17337 :                                      qs_kind_set=qs_kind_set, ires=ires)
     151              : 
     152        17337 :                CALL cp_print_key_finished_output(ires, logger, dft_section, TRIM(keys(ikey)))
     153              :             END IF
     154              :          END DO
     155              :       END IF
     156              : 
     157       161937 :       CALL timestop(handle)
     158              : 
     159       161937 :    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         9681 :    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         9681 :       CALL timeset(routineN, handle)
     185         9681 :       logger => cp_get_default_logger()
     186         9681 :       IF (logger%para_env%is_source()) THEN
     187         4968 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     188              :       ELSE
     189              :          unit_nr = -1
     190              :       END IF
     191              : 
     192         9681 :       project_name = logger%iter_info%project_name
     193         9681 :       CALL section_vals_val_get(dft_section, "SCF%PRINT%DM_RESTART_WRITE", l_val=do_dm_restart)
     194         9681 :       NULLIFY (matrix_p_tmp)
     195              : 
     196         9681 :       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         9681 :       CALL timestop(handle)
     220              : 
     221         9681 :    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        17453 :    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        17453 :       INTEGER, DIMENSION(:), POINTER                     :: nset_info, nshell
     300        17453 :       INTEGER, DIMENSION(:, :), POINTER                  :: l, nshell_info
     301        17453 :       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        17453 :       CALL timeset(routineN, handle)
     306        17453 :       nspin = SIZE(mo_array)
     307        17453 :       nao = mo_array(1)%nao
     308              : 
     309        17453 :       IF (ires > 0) THEN
     310              :          !     *** create some info about the basis set first ***
     311         8898 :          natom = SIZE(particle_set, 1)
     312         8898 :          nset_max = 0
     313         8898 :          nshell_max = 0
     314              : 
     315        56525 :          DO iatom = 1, natom
     316        47627 :             NULLIFY (orb_basis_set, dftb_parameter)
     317        47627 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     318              :             CALL get_qs_kind(qs_kind_set(ikind), &
     319        47627 :                              basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
     320       104152 :             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        40089 :                                       l=l)
     325        40089 :                nset_max = MAX(nset_max, nset)
     326       113991 :                DO iset = 1, nset
     327       113991 :                   nshell_max = MAX(nshell_max, nshell(iset))
     328              :                END DO
     329         7538 :             ELSEIF (ASSOCIATED(dftb_parameter)) THEN
     330         7537 :                CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
     331         7537 :                nset_max = MAX(nset_max, 1)
     332         7537 :                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        44490 :          ALLOCATE (nso_info(nshell_max, nset_max, natom))
     340       326547 :          nso_info(:, :, :) = 0
     341              : 
     342        35592 :          ALLOCATE (nshell_info(nset_max, natom))
     343       148076 :          nshell_info(:, :) = 0
     344              : 
     345        26694 :          ALLOCATE (nset_info(natom))
     346        56525 :          nset_info(:) = 0
     347              : 
     348        56525 :          DO iatom = 1, natom
     349        47627 :             NULLIFY (orb_basis_set, dftb_parameter)
     350        47627 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     351              :             CALL get_qs_kind(qs_kind_set(ikind), &
     352        47627 :                              basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
     353       104152 :             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        40089 :                                       l=l)
     358        40089 :                nset_info(iatom) = nset
     359       113991 :                DO iset = 1, nset
     360        73902 :                   nshell_info(iset, iatom) = nshell(iset)
     361       217562 :                   DO ishell = 1, nshell(iset)
     362       103571 :                      lshell = l(ishell, iset)
     363       177473 :                      nso_info(ishell, iset, iatom) = nso(lshell)
     364              :                   END DO
     365              :                END DO
     366         7538 :             ELSEIF (ASSOCIATED(dftb_parameter)) THEN
     367         7537 :                CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
     368         7537 :                nset_info(iatom) = 1
     369         7537 :                nshell_info(1, iatom) = lmax + 1
     370        17656 :                DO ishell = 1, lmax + 1
     371        10119 :                   lshell = ishell - 1
     372        17656 :                   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         8898 :          WRITE (ires) natom, nspin, nao, nset_max, nshell_max
     381        56525 :          WRITE (ires) nset_info
     382       148076 :          WRITE (ires) nshell_info
     383       326547 :          WRITE (ires) nso_info
     384              : 
     385         8898 :          DEALLOCATE (nset_info)
     386              : 
     387         8898 :          DEALLOCATE (nshell_info)
     388              : 
     389         8898 :          DEALLOCATE (nso_info)
     390              :       END IF
     391              : 
     392              :       ! use the scalapack block size as a default for buffering columns
     393        17453 :       CALL cp_fm_get_info(mo_array(1)%mo_coeff, ncol_block=max_block)
     394        37780 :       DO ispin = 1, nspin
     395        20327 :          nmo = mo_array(ispin)%nmo
     396        20327 :          IF ((ires > 0) .AND. (nmo > 0)) THEN
     397        10230 :             WRITE (ires) nmo, &
     398        10230 :                mo_array(ispin)%homo, &
     399        10230 :                mo_array(ispin)%lfomo, &
     400        20460 :                mo_array(ispin)%nelectron
     401        99233 :             WRITE (ires) mo_array(ispin)%eigenvalues(1:nmo), &
     402       109463 :                mo_array(ispin)%occupation_numbers(1:nmo)
     403              :          END IF
     404        37780 :          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        20189 :             CALL cp_fm_write_unformatted(mo_array(ispin)%mo_coeff, ires)
     410              :          END IF
     411              :       END DO
     412              : 
     413        17453 :       CALL timestop(handle)
     414              : 
     415        17453 :    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         1204 :    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          602 :       my_kp = .FALSE.
     439          602 :       my_xas = .FALSE.
     440          602 :       my_rtp = .FALSE.
     441          602 :       IF (PRESENT(kp)) my_kp = kp
     442          602 :       IF (PRESENT(xas)) my_xas = xas
     443          602 :       IF (PRESENT(rtp)) my_rtp = rtp
     444              : 
     445          602 :       exist = .FALSE.
     446          602 :       CALL section_vals_val_get(section, "WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
     447          602 :       IF (n_rep_val > 0) THEN
     448          462 :          CALL section_vals_val_get(section, "WFN_RESTART_FILE_NAME", c_val=filename)
     449              :       ELSE
     450          140 :          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          136 :          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          133 :          ELSE IF (my_kp) THEN
     461              :             ! try to read from the filename that is generated automatically from the printkey
     462           12 :             print_key => section_vals_get_subs_vals(section, "SCF%PRINT%RESTART")
     463              :             filename = cp_print_key_generate_filename(logger, print_key, &
     464           12 :                                                       extension=".kp", my_local=.FALSE.)
     465              :          ELSE
     466              :             ! try to read from the filename that is generated automatically from the printkey
     467          121 :             print_key => section_vals_get_subs_vals(section, "SCF%PRINT%RESTART")
     468              :             filename = cp_print_key_generate_filename(logger, print_key, &
     469          121 :                                                       extension=".wfn", my_local=.FALSE.)
     470              :          END IF
     471              :       END IF
     472          602 :       IF (.NOT. my_xas) THEN
     473          596 :          INQUIRE (FILE=filename, exist=exist)
     474              :       END IF
     475              : 
     476          602 :    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          521 :    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          521 :       CALL timeset(routineN, handle)
     516          521 :       logger => cp_get_default_logger()
     517          521 :       my_cdft = .FALSE.
     518          521 :       IF (PRESENT(cdft)) my_cdft = cdft
     519          521 :       my_out_unit = -1
     520          521 :       IF (PRESENT(out_unit)) my_out_unit = out_unit
     521              : 
     522          521 :       nspin = SIZE(mo_array)
     523          521 :       restart_unit = -1
     524              : 
     525          521 :       IF (para_env%is_source()) THEN
     526              : 
     527          279 :          natom = SIZE(particle_set, 1)
     528          279 :          CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
     529          279 :          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          279 :                         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          521 :                                 rst_unit=restart_unit, multiplicity=multiplicity, natom_mismatch=natom_mismatch)
     545              : 
     546          521 :       IF (PRESENT(natom_mismatch)) THEN
     547              :          ! read_mos_restart_low only the io_node returns natom_mismatch, must broadcast it
     548          491 :          CALL para_env%bcast(natom_mismatch)
     549          491 :          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          521 :       IF (para_env%is_source()) THEN
     558          279 :          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          279 :          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          521 :       IF (.NOT. my_cdft) THEN
     567         1396 :          DO ispin = 1, nspin
     568              :             CALL write_mo_set_to_output_unit(mo_array(ispin), atomic_kind_set, qs_kind_set, &
     569         1396 :                                              particle_set, dft_section, 4, 0, final_mos=.FALSE.)
     570              :          END DO
     571              :       END IF
     572              : 
     573          521 :       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          585 :    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          585 :       INTEGER, DIMENSION(:), POINTER                     :: nset_info, nshell
     685          585 :       INTEGER, DIMENSION(:, :), POINTER                  :: l, nshell_info
     686          585 :       INTEGER, DIMENSION(:, :, :), POINTER               :: nso_info, offset_info
     687              :       LOGICAL                                            :: minbas, natom_match, use_this
     688          585 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eig_read, occ_read
     689          585 :       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         1170 :       logger => cp_get_default_logger()
     695              : 
     696          585 :       nspin = SIZE(mos)
     697          585 :       nao = mos(1)%nao
     698          585 :       my_mult = 0
     699          585 :       IF (PRESENT(multiplicity)) my_mult = multiplicity
     700              : 
     701          585 :       IF (para_env%is_source()) THEN
     702          311 :          READ (rst_unit) natom_read, nspin_read, nao_read, nset_max, nshell_max
     703          311 :          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          307 :             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          307 :             IF (nspin_read > nspin) THEN
     715            0 :                CPABORT("Reducing nspin is not possible. ")
     716              :             END IF
     717              :          END IF
     718              : 
     719          311 :          natom_match = (natom_read == natom)
     720              : 
     721          311 :          IF (natom_match) THEN ! actually do the read read
     722              : 
     723              :             ! Let's make it possible to change the basis set
     724         1555 :             ALLOCATE (nso_info(nshell_max, nset_max, natom_read))
     725         1244 :             ALLOCATE (nshell_info(nset_max, natom_read))
     726          933 :             ALLOCATE (nset_info(natom_read))
     727         1244 :             ALLOCATE (offset_info(nshell_max, nset_max, natom_read))
     728              : 
     729          311 :             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         1140 :             READ (rst_unit) nset_info
     737         2865 :             READ (rst_unit) nshell_info
     738         6943 :             READ (rst_unit) nso_info
     739              : 
     740          311 :             i = 1
     741         1140 :             DO iatom = 1, natom
     742         2681 :                DO iset = 1, nset_info(iatom)
     743         4900 :                   DO ishell = 1, nshell_info(iset, iatom)
     744         2530 :                      offset_info(ishell, iset, iatom) = i
     745         4071 :                      i = i + nso_info(ishell, iset, iatom)
     746              :                   END DO
     747              :                END DO
     748              :             END DO
     749              : 
     750          933 :             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          585 :       CALL para_env%bcast(natom_match)
     757          585 :       IF (PRESENT(natom_mismatch)) natom_mismatch = .NOT. natom_match
     758              :       ! handle natom_match false
     759          585 :       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          585 :       CALL para_env%bcast(nspin_read)
     770              : 
     771         1755 :       ALLOCATE (vecbuffer(1, nao))
     772              : 
     773         1602 :       DO ispin = 1, nspin
     774              : 
     775         1017 :          nmo = mos(ispin)%nmo
     776         1017 :          homo = mos(ispin)%homo
     777         5340 :          mos(ispin)%eigenvalues(:) = 0.0_dp
     778         5340 :          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         2749 :             eig_read = 0.0_dp
     785         2749 :             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         2725 :             mos(ispin)%eigenvalues(1:nmo) = eig_read(1:nmo)
     799         2725 :             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         9663 :          CALL para_env%bcast(mos(ispin)%eigenvalues)
     823         9663 :          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         5288 :             DO i = 1, nmo
     840         4279 :                IF (para_env%is_source()) THEN
     841       146193 :                   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         2185 :                   irow = 1
     845         9774 :                   DO iatom = 1, natom
     846         7589 :                      NULLIFY (orb_basis_set, dftb_parameter, l, nshell)
     847         7589 :                      CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     848              :                      CALL get_qs_kind(qs_kind_set(ikind), &
     849         7589 :                                       basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
     850         7589 :                      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         7541 :                                                l=l)
     855         7541 :                         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         7589 :                      use_this = .TRUE.
     867         7589 :                      iset_read = 1
     868        35528 :                      DO iset = 1, nset
     869        18165 :                         ishell_read = 1
     870        18165 :                         IF (minbas) THEN
     871           48 :                            nnshell = lmax + 1
     872              :                         ELSE
     873        18117 :                            nnshell = nshell(iset)
     874              :                         END IF
     875        57458 :                         DO ishell = 1, nnshell
     876        31704 :                            IF (minbas) THEN
     877           72 :                               lshell = ishell - 1
     878              :                            ELSE
     879        31632 :                               lshell = l(ishell, iset)
     880              :                            END IF
     881        31704 :                            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        31680 :                               IF (nso(lshell) == nso_info(ishell_read, iset_read, iatom)) THEN
     884        31680 :                                  offset_read = offset_info(ishell_read, iset_read, iatom)
     885        31680 :                                  ishell_read = ishell_read + 1
     886        31680 :                                  IF (ishell_read > nshell_info(iset, iatom)) THEN
     887        18161 :                                     ishell_read = 1
     888        18161 :                                     iset_read = iset_read + 1
     889              :                                  END IF
     890              :                               ELSE
     891              :                                  use_this = .FALSE.
     892              :                               END IF
     893              :                            END IF
     894       103852 :                            DO iso = 1, nso(lshell)
     895        72148 :                               IF (use_this) THEN
     896        72004 :                                  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        72004 :                                     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       103852 :                               irow = irow + 1
     905              :                            END DO
     906        49869 :                            use_this = .TRUE.
     907              :                         END DO
     908              :                      END DO
     909              :                   END DO
     910              : 
     911              :                ELSE
     912              : 
     913        72764 :                   vecbuffer(1, :) = 0.0_dp
     914              : 
     915              :                END IF
     916              : 
     917       575551 :                CALL para_env%bcast(vecbuffer)
     918              :                CALL cp_fm_set_submatrix(mos(ispin)%mo_coeff, &
     919         5288 :                                         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         1602 :          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          585 :       DEALLOCATE (vecbuffer)
     955              : 
     956          585 :       IF (para_env%is_source()) THEN
     957          311 :          DEALLOCATE (vecbuffer_read)
     958          311 :          DEALLOCATE (offset_info)
     959          311 :          DEALLOCATE (nso_info)
     960          311 :          DEALLOCATE (nshell_info)
     961          311 :          DEALLOCATE (nset_info)
     962              :       END IF
     963              : 
     964         1170 :    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         8985 :    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         8985 :       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         8985 :       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         8985 :       INTEGER, DIMENSION(:), POINTER                     :: mo_index_range, nshell
    1027         8985 :       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         8985 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: mo_eigenvalues, mo_occupation_numbers
    1034         8985 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cmatrix, smatrix
    1035         8985 :       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         8985 :       NULLIFY (bcgf_symbol)
    1042         8985 :       NULLIFY (bsgf_symbol)
    1043         8985 :       NULLIFY (logger)
    1044         8985 :       NULLIFY (mo_index_range)
    1045         8985 :       NULLIFY (nshell)
    1046         8985 :       NULLIFY (mo_coeff)
    1047              : 
    1048        17970 :       logger => cp_get_default_logger()
    1049         8985 :       ionode = logger%para_env%is_source()
    1050         8985 :       CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVALUES", l_val=print_eigvals)
    1051         8985 :       CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVECTORS", l_val=print_eigvecs)
    1052         8985 :       CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occup)
    1053         8985 :       CALL section_vals_val_get(dft_section, "PRINT%MO%CARTESIAN", l_val=print_cartesian)
    1054         8985 :       CALL section_vals_val_get(dft_section, "PRINT%MO%MO_INDEX_RANGE", i_vals=mo_index_range)
    1055         8985 :       CALL section_vals_val_get(dft_section, "PRINT%MO%NDIGITS", i_val=after)
    1056         8985 :       after = MIN(MAX(after, 1), 16)
    1057              : 
    1058              :       ! Do we print the final MO information after SCF convergence is reached (default: no)
    1059         8985 :       IF (PRESENT(final_mos)) THEN
    1060         8977 :          my_final = final_mos
    1061              :       ELSE
    1062              :          my_final = .FALSE.
    1063              :       END IF
    1064              : 
    1065              :       ! complex MOS for RTP, no eigenvalues
    1066         8985 :       my_rtp = .FALSE.
    1067         8985 :       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         9898 :                                                           "PRINT%MO"), cp_p_file) .OR. my_final
    1076              :       END IF
    1077              : 
    1078         8985 :       IF ((.NOT. should_output) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup))) RETURN
    1079              : 
    1080         8062 :       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         8054 :          scf_step = MAX(0, logger%iter_info%iteration(logger%iter_info%n_rlevel) - 1)
    1092              :       END IF
    1093              : 
    1094         8062 :       IF (.NOT. my_final) THEN
    1095         6928 :          IF (.NOT. my_rtp) THEN
    1096         6920 :             step_string = " AFTER SCF STEP"
    1097              :          ELSE
    1098            8 :             step_string = " AFTER RTP STEP"
    1099              :          END IF
    1100              :       END IF
    1101              : 
    1102         8062 :       IF (PRESENT(solver_method)) THEN
    1103         7908 :          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         8062 :                       nmo=nmo)
    1118         8062 :       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         8042 :          numo = 0
    1125              :       END IF
    1126        24136 :       ALLOCATE (mo_eigenvalues(nmo))
    1127        58302 :       mo_eigenvalues(:) = 0.0_dp
    1128        58192 :       mo_eigenvalues(1:nmo - numo) = eigenvalues(1:nmo - numo)
    1129        16074 :       ALLOCATE (mo_occupation_numbers(nmo))
    1130        58302 :       mo_occupation_numbers(:) = 0.0_dp
    1131        58192 :       mo_occupation_numbers(1:nmo - numo) = occupation_numbers(1:nmo - numo)
    1132         8062 :       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         8062 :       IF (print_eigvecs) THEN
    1139        25814 :          ALLOCATE (smatrix(nao, nmo))
    1140         6466 :          CALL cp_fm_get_submatrix(mo_coeff, smatrix(1:nao, 1:nmo - numo))
    1141         6466 :          IF (numo > 0) THEN
    1142           14 :             CALL cp_fm_get_submatrix(umo_coeff, smatrix(1:nao, nmo - numo + 1:nmo))
    1143              :          END IF
    1144         6466 :          IF (.NOT. ionode) THEN
    1145         3233 :             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         8062 :                                 extension=".MOLog")
    1152              : 
    1153         8062 :       IF (iw > 0) THEN
    1154              : 
    1155         4031 :          CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
    1156         4031 :          CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
    1157              : 
    1158              :          ! Definition of the variable formats
    1159              : 
    1160         4031 :          fmtstr1 = "(T2,A,21X,  (  X,I5,  X))"
    1161         4031 :          fmtstr2 = "(T2,A,21X,  (1X,F  .  ))"
    1162         4031 :          fmtstr3 = "(T2,A,I5,1X,I5,1X,A,1X,A6,  (1X,F  .  ))"
    1163              : 
    1164         4031 :          width = before + after + 3
    1165         4031 :          ncol = INT(56/width)
    1166              : 
    1167         4031 :          right = MAX((after - 2), 1)
    1168         4031 :          left = width - right - 5
    1169              : 
    1170         4031 :          WRITE (UNIT=fmtstr1(11:12), FMT="(I2)") ncol
    1171         4031 :          WRITE (UNIT=fmtstr1(14:15), FMT="(I2)") left
    1172         4031 :          WRITE (UNIT=fmtstr1(21:22), FMT="(I2)") right
    1173              : 
    1174         4031 :          WRITE (UNIT=fmtstr2(11:12), FMT="(I2)") ncol
    1175         4031 :          WRITE (UNIT=fmtstr2(18:19), FMT="(I2)") width - 1
    1176         4031 :          WRITE (UNIT=fmtstr2(21:22), FMT="(I2)") after
    1177              : 
    1178         4031 :          WRITE (UNIT=fmtstr3(27:28), FMT="(I2)") ncol
    1179         4031 :          WRITE (UNIT=fmtstr3(34:35), FMT="(I2)") width - 1
    1180         4031 :          WRITE (UNIT=fmtstr3(37:38), FMT="(I2)") after
    1181              : 
    1182         4031 :          IF (my_final .OR. (my_solver_method == "TD")) THEN
    1183         2539 :             energy_str = "EIGENVALUES"
    1184         2539 :             vector_str = "EIGENVECTORS"
    1185              :          ELSE
    1186         1492 :             energy_str = "ENERGIES"
    1187         1492 :             vector_str = "COEFFICIENTS"
    1188              :          END IF
    1189              : 
    1190         4031 :          IF (my_rtp) THEN
    1191            4 :             energy_str = "ZEROS"
    1192            4 :             vector_str = TRIM(reim)//" RTP COEFFICIENTS"
    1193              :          END IF
    1194              : 
    1195         4031 :          IF (print_eigvecs) THEN
    1196              : 
    1197         3233 :             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         3148 :                orbital_str = "SPHERICAL"
    1249              : 
    1250              :             END IF ! print_cartesian
    1251              : 
    1252              :             name = TRIM(energy_str)//", OCCUPATION NUMBERS, AND "// &
    1253         3233 :                    TRIM(orbital_str)//" "//TRIM(vector_str)
    1254              : 
    1255         3233 :             IF (.NOT. my_final) &
    1256         2889 :                WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
    1257              : 
    1258          798 :          ELSE IF (print_occup .OR. print_eigvals) THEN
    1259          798 :             name = TRIM(energy_str)//" AND OCCUPATION NUMBERS"
    1260              : 
    1261          798 :             IF (.NOT. my_final) &
    1262          575 :                WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
    1263              :          END IF ! print_eigvecs
    1264              : 
    1265              :          ! Print headline
    1266         4031 :          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         3707 :          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         3606 :                "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         4031 :              (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         3860 :             first_mo = 1
    1288         3860 :             last_mo = nmo
    1289              :          END IF
    1290              : 
    1291         4031 :          IF (print_eigvecs) THEN
    1292              : 
    1293              :             ! Print full MO information
    1294              : 
    1295         6912 :             DO icol = first_mo, last_mo, ncol
    1296              : 
    1297         3679 :                from = icol
    1298         3679 :                to = MIN((from + ncol - 1), last_mo)
    1299              : 
    1300         3679 :                WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1301              :                WRITE (UNIT=iw, FMT=fmtstr1) &
    1302        17555 :                   "MO|", (jcol, jcol=from, to)
    1303              :                WRITE (UNIT=iw, FMT=fmtstr2) &
    1304         3679 :                   "MO|", (mo_eigenvalues(jcol), jcol=from, to)
    1305         3679 :                WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1306              :                WRITE (UNIT=iw, FMT=fmtstr2) &
    1307         3679 :                   "MO|", (mo_occupation_numbers(jcol), jcol=from, to)
    1308         3679 :                WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1309              : 
    1310         3679 :                irow = 1
    1311              : 
    1312        14619 :                DO iatom = 1, natom
    1313              : 
    1314         7707 :                   IF (iatom /= 1) WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1315              : 
    1316         7707 :                   NULLIFY (orb_basis_set, dftb_parameter)
    1317              :                   CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
    1318         7707 :                                        element_symbol=element_symbol, kind_number=ikind)
    1319              :                   CALL get_qs_kind(qs_kind_set(ikind), &
    1320              :                                    basis_set=orb_basis_set, &
    1321         7707 :                                    dftb_parameter=dftb_parameter)
    1322              : 
    1323        19093 :                   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         6855 :                      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         6855 :                                                sgf_symbol=bsgf_symbol)
    1373         6855 :                         isgf = 1
    1374        20015 :                         DO iset = 1, nset
    1375        35573 :                            DO ishell = 1, nshell(iset)
    1376        15558 :                               lshell = l(ishell, iset)
    1377        66194 :                               DO iso = 1, nso(lshell)
    1378              :                                  WRITE (UNIT=iw, FMT=fmtstr3) &
    1379        37476 :                                     "MO|", irow, iatom, ADJUSTR(element_symbol), bsgf_symbol(isgf), &
    1380        74952 :                                     (smatrix(irow, jcol), jcol=from, to)
    1381        37476 :                                  isgf = isgf + 1
    1382        53034 :                                  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         3233 :             WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1413              : 
    1414              :             ! Release work storage
    1415              : 
    1416         3233 :             DEALLOCATE (smatrix)
    1417         3233 :             IF (print_cartesian) THEN
    1418           85 :                DEALLOCATE (cmatrix)
    1419              :             END IF
    1420              : 
    1421          798 :          ELSE IF (print_occup .OR. print_eigvals) THEN
    1422              : 
    1423          798 :             WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
    1424          798 :             fmtstr4 = "(T2,A,I7,3(1X,F22.  ))"
    1425          798 :             WRITE (UNIT=fmtstr4(19:20), FMT="(I2)") after
    1426          798 :             IF (my_final .OR. (my_solver_method == "TD")) THEN
    1427              :                WRITE (UNIT=iw, FMT="(A)") &
    1428          798 :                   " 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        11998 :             DO imo = first_mo, last_mo
    1434              :                WRITE (UNIT=iw, FMT=fmtstr4) &
    1435        11200 :                   "MO|", imo, mo_eigenvalues(imo), &
    1436        11200 :                   mo_eigenvalues(imo)*evolt, &
    1437        23198 :                   mo_occupation_numbers(imo)
    1438              :             END DO
    1439          798 :             fmtstr5 = "(A,T59,F22.  )"
    1440          798 :             WRITE (UNIT=fmtstr5(12:13), FMT="(I2)") after
    1441              :             WRITE (UNIT=iw, FMT=fmtstr5) &
    1442          798 :                " MO| Sum:", accurate_sum(mo_occupation_numbers(:))
    1443              : 
    1444              :          END IF ! print_eigvecs
    1445              : 
    1446         4031 :          IF (.NOT. my_rtp) THEN
    1447         4027 :             fmtstr6 = "(A,T18,F17.  ,A,T41,F17.  ,A)"
    1448         4027 :             WRITE (UNIT=fmtstr6(12:13), FMT="(I2)") after
    1449         4027 :             WRITE (UNIT=fmtstr6(25:26), FMT="(I2)") after
    1450              :             WRITE (UNIT=iw, FMT=fmtstr6) &
    1451         4027 :                " MO| E(Fermi):", mo_set%mu, " a.u.", mo_set%mu*evolt, " eV"
    1452              :          END IF
    1453         4031 :          IF ((homo > 0) .AND. .NOT. my_rtp) THEN
    1454         4002 :             IF ((mo_occupation_numbers(homo) == maxocc) .AND. (last_mo > homo)) THEN
    1455              :                gap = mo_eigenvalues(homo + 1) - &
    1456          351 :                      mo_eigenvalues(homo)
    1457              :                WRITE (UNIT=iw, FMT=fmtstr6) &
    1458          351 :                   " MO| Band gap:", gap, " a.u.", gap*evolt, " eV"
    1459              :             END IF
    1460              :          END IF
    1461         4031 :          WRITE (UNIT=iw, FMT="(A)") ""
    1462              : 
    1463              :       END IF ! iw
    1464              : 
    1465         8062 :       IF (ALLOCATED(mo_eigenvalues)) DEALLOCATE (mo_eigenvalues)
    1466         8062 :       IF (ALLOCATED(mo_occupation_numbers)) DEALLOCATE (mo_occupation_numbers)
    1467              : 
    1468              :       CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%MO", &
    1469         8062 :                                         ignore_should_output=should_output)
    1470              : 
    1471        35017 :    END SUBROUTINE write_mo_set_to_output_unit
    1472              : 
    1473              : END MODULE qs_mo_io
        

Generated by: LCOV version 2.0-1