LCOV - code coverage report
Current view: top level - src - qs_pdos.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 65.9 % 464 306
Test Date: 2025-07-25 12:55:17 Functions: 33.3 % 6 2

            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 Calculation and writing of projected density of  states
      10              : !>         The DOS is computed per angular momentum and per kind
      11              : !> \par History
      12              : !>      -
      13              : !> \author Marcella (29.02.2008,MK)
      14              : ! **************************************************************************************************
      15              : MODULE qs_pdos
      16              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17              :                                               get_atomic_kind,&
      18              :                                               get_atomic_kind_set
      19              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      20              :                                               gto_basis_set_type
      21              :    USE cell_types,                      ONLY: cell_type,&
      22              :                                               pbc
      23              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      24              :    USE cp_control_types,                ONLY: dft_control_type
      25              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      26              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm
      27              :    USE cp_fm_diag,                      ONLY: cp_fm_power
      28              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      29              :                                               cp_fm_struct_release,&
      30              :                                               cp_fm_struct_type
      31              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      32              :                                               cp_fm_get_info,&
      33              :                                               cp_fm_get_submatrix,&
      34              :                                               cp_fm_init_random,&
      35              :                                               cp_fm_release,&
      36              :                                               cp_fm_type
      37              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      38              :                                               cp_logger_get_default_io_unit,&
      39              :                                               cp_logger_type,&
      40              :                                               cp_to_string
      41              :    USE cp_output_handling,              ONLY: cp_p_file,&
      42              :                                               cp_print_key_finished_output,&
      43              :                                               cp_print_key_should_output,&
      44              :                                               cp_print_key_unit_nr
      45              :    USE input_section_types,             ONLY: section_vals_get,&
      46              :                                               section_vals_get_subs_vals,&
      47              :                                               section_vals_type,&
      48              :                                               section_vals_val_get
      49              :    USE kinds,                           ONLY: default_string_length,&
      50              :                                               dp
      51              :    USE memory_utilities,                ONLY: reallocate
      52              :    USE message_passing,                 ONLY: mp_para_env_type
      53              :    USE orbital_pointers,                ONLY: nso,&
      54              :                                               nsoset
      55              :    USE orbital_symbols,                 ONLY: l_sym,&
      56              :                                               sgf_symbol
      57              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      58              :    USE particle_types,                  ONLY: particle_type
      59              :    USE preconditioner_types,            ONLY: preconditioner_type
      60              :    USE pw_env_types,                    ONLY: pw_env_get,&
      61              :                                               pw_env_type
      62              :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      63              :                                               pw_pool_type
      64              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      65              :                                               pw_r3d_rs_type
      66              :    USE qs_collocate_density,            ONLY: calculate_wavefunction
      67              :    USE qs_environment_types,            ONLY: get_qs_env,&
      68              :                                               qs_environment_type
      69              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      70              :                                               get_qs_kind_set,&
      71              :                                               qs_kind_type
      72              :    USE qs_mo_methods,                   ONLY: calculate_subspace_eigenvalues
      73              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      74              :                                               mo_set_type
      75              :    USE qs_ot_eigensolver,               ONLY: ot_eigensolver
      76              :    USE scf_control_types,               ONLY: scf_control_type
      77              : #include "./base/base_uses.f90"
      78              : 
      79              :    IMPLICIT NONE
      80              : 
      81              :    PRIVATE
      82              : 
      83              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_pdos'
      84              : 
      85              : ! **************************************************************************************************
      86              :    ! *** Public subroutines ***
      87              : 
      88              :    PUBLIC :: calculate_projected_dos
      89              : 
      90              :    TYPE ldos_type
      91              :       INTEGER :: maxl = -1, nlist = -1
      92              :       LOGICAL :: separate_components = .FALSE.
      93              :       INTEGER, DIMENSION(:), POINTER           :: list_index => NULL()
      94              :       REAL(KIND=dp), DIMENSION(:, :), &
      95              :          POINTER                                :: pdos_array => NULL()
      96              :    END TYPE ldos_type
      97              : 
      98              :    TYPE r_ldos_type
      99              :       INTEGER :: nlist = -1, npoints = -1
     100              :       INTEGER, DIMENSION(:, :), POINTER :: index_grid_local => NULL()
     101              :       INTEGER, DIMENSION(:), POINTER           :: list_index => NULL()
     102              :       REAL(KIND=dp), DIMENSION(:), POINTER     :: x_range => NULL(), y_range => NULL(), z_range => NULL()
     103              :       REAL(KIND=dp), DIMENSION(:), POINTER     :: eval_range => NULL()
     104              :       REAL(KIND=dp), DIMENSION(:), &
     105              :          POINTER                                 :: pdos_array => NULL()
     106              :    END TYPE r_ldos_type
     107              : 
     108              :    TYPE ldos_p_type
     109              :       TYPE(ldos_type), POINTER :: ldos => NULL()
     110              :    END TYPE ldos_p_type
     111              : 
     112              :    TYPE r_ldos_p_type
     113              :       TYPE(r_ldos_type), POINTER :: ldos => NULL()
     114              :    END TYPE r_ldos_p_type
     115              : CONTAINS
     116              : 
     117              : ! **************************************************************************************************
     118              : !> \brief   Compute and write projected density of states
     119              : !> \param mo_set ...
     120              : !> \param atomic_kind_set ...
     121              : !> \param qs_kind_set ...
     122              : !> \param particle_set ...
     123              : !> \param qs_env ...
     124              : !> \param dft_section ...
     125              : !> \param ispin ...
     126              : !> \param xas_mittle ...
     127              : !> \param external_matrix_shalf ...
     128              : !> \date    26.02.2008
     129              : !> \par History:
     130              : !>       - Added optional external matrix_shalf to avoid recomputing it (A. Bussy, 09.2019)
     131              : !> \par Variables
     132              : !>       -
     133              : !>       -
     134              : !> \author  MI
     135              : !> \version 1.0
     136              : ! **************************************************************************************************
     137           60 :    SUBROUTINE calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, &
     138              :                                       dft_section, ispin, xas_mittle, external_matrix_shalf)
     139              : 
     140              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     141              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     142              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     143              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     144              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     145              :       TYPE(section_vals_type), POINTER                   :: dft_section
     146              :       INTEGER, INTENT(IN), OPTIONAL                      :: ispin
     147              :       CHARACTER(LEN=default_string_length), INTENT(IN), &
     148              :          OPTIONAL                                        :: xas_mittle
     149              :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL, TARGET     :: external_matrix_shalf
     150              : 
     151              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos'
     152              : 
     153              :       CHARACTER(LEN=16)                                  :: fmtstr2
     154              :       CHARACTER(LEN=27)                                  :: fmtstr1
     155           60 :       CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:, :, :)  :: tmp_str
     156              :       CHARACTER(LEN=default_string_length)               :: kind_name, my_act, my_mittle, my_pos, &
     157              :                                                             spin(2)
     158              :       CHARACTER(LEN=default_string_length), &
     159           60 :          ALLOCATABLE, DIMENSION(:)                       :: ldos_index, r_ldos_index
     160              :       INTEGER :: handle, homo, i, iatom, ikind, il, ildos, im, imo, in_x, in_y, in_z, ir, irow, &
     161              :          iset, isgf, ishell, iso, iterstep, iw, j, jx, jy, jz, k, lcomponent, lshell, maxl, &
     162              :          maxlgto, my_spin, n_dependent, n_r_ldos, n_rep, nao, natom, ncol_global, nkind, nldos, &
     163              :          nlumo, nmo, np_tot, npoints, nrow_global, nset, nsgf, nvirt, out_each, output_unit
     164           60 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: firstrow
     165           60 :       INTEGER, DIMENSION(:), POINTER                     :: list, nshell
     166           60 :       INTEGER, DIMENSION(:, :), POINTER                  :: bo, l
     167              :       LOGICAL                                            :: append, calc_matsh, do_ldos, do_r_ldos, &
     168              :                                                             do_virt, ionode, separate_components, &
     169              :                                                             should_output
     170           60 :       LOGICAL, DIMENSION(:, :), POINTER                  :: read_r
     171              :       REAL(KIND=dp)                                      :: dh(3, 3), dvol, e_fermi, r(3), r_vec(3), &
     172              :                                                             ratom(3)
     173           60 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eigenvalues, evals_virt, &
     174           60 :                                                             occupation_numbers
     175           60 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
     176           60 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: pdos_array
     177              :       TYPE(cell_type), POINTER                           :: cell
     178              :       TYPE(cp_blacs_env_type), POINTER                   :: context
     179              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     180              :       TYPE(cp_fm_type)                                   :: matrix_shalfc, matrix_work, mo_virt
     181              :       TYPE(cp_fm_type), POINTER                          :: matrix_shalf, mo_coeff
     182              :       TYPE(cp_logger_type), POINTER                      :: logger
     183           60 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: s_matrix
     184              :       TYPE(dft_control_type), POINTER                    :: dft_control
     185              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     186           60 :       TYPE(ldos_p_type), DIMENSION(:), POINTER           :: ldos_p
     187              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     188              :       TYPE(pw_c1d_gs_type)                               :: wf_g
     189              :       TYPE(pw_env_type), POINTER                         :: pw_env
     190           60 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     191              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     192              :       TYPE(pw_r3d_rs_type)                               :: wf_r
     193           60 :       TYPE(r_ldos_p_type), DIMENSION(:), POINTER         :: r_ldos_p
     194              :       TYPE(section_vals_type), POINTER                   :: ldos_section
     195              : 
     196           60 :       NULLIFY (logger)
     197          120 :       logger => cp_get_default_logger()
     198              :       ionode = logger%para_env%is_source()
     199              :       should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
     200           60 :                                                        "PRINT%PDOS"), cp_p_file)
     201           60 :       output_unit = cp_logger_get_default_io_unit(logger)
     202              : 
     203           60 :       spin(1) = "ALPHA"
     204           60 :       spin(2) = "BETA"
     205           60 :       IF ((.NOT. should_output)) RETURN
     206              : 
     207           60 :       NULLIFY (context, s_matrix, orb_basis_set, para_env, pdos_array)
     208           60 :       NULLIFY (eigenvalues, fm_struct_tmp, mo_coeff, vecbuffer)
     209           60 :       NULLIFY (ldos_section, list, cell, pw_env, auxbas_pw_pool, evals_virt)
     210           60 :       NULLIFY (occupation_numbers, ldos_p, r_ldos_p, dft_control, occupation_numbers)
     211              : 
     212           60 :       CALL timeset(routineN, handle)
     213           60 :       iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
     214              : 
     215           60 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
     216           30 :          " Calculate PDOS at iteration step ", iterstep
     217              :       CALL get_qs_env(qs_env=qs_env, &
     218           60 :                       matrix_s=s_matrix)
     219              : 
     220           60 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
     221           60 :       CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf, maxlgto=maxlgto)
     222           60 :       nkind = SIZE(atomic_kind_set)
     223              : 
     224              :       CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo, &
     225           60 :                       mu=e_fermi)
     226              :       CALL cp_fm_get_info(mo_coeff, &
     227              :                           context=context, para_env=para_env, &
     228              :                           nrow_global=nrow_global, &
     229           60 :                           ncol_global=ncol_global)
     230              : 
     231           60 :       CALL section_vals_val_get(dft_section, "PRINT%PDOS%OUT_EACH_MO", i_val=out_each)
     232           60 :       IF (out_each == -1) out_each = nao + 1
     233           60 :       CALL section_vals_val_get(dft_section, "PRINT%PDOS%NLUMO", i_val=nlumo)
     234           60 :       IF (nlumo == -1) nlumo = nao - homo
     235           60 :       do_virt = (nlumo > (nmo - homo))
     236           60 :       nvirt = nlumo - (nmo - homo)
     237              :       ! Generate virtual orbitals
     238           60 :       IF (do_virt) THEN
     239           14 :          IF (PRESENT(ispin)) THEN
     240            8 :             my_spin = ispin
     241              :          ELSE
     242            6 :             my_spin = 1
     243              :          END IF
     244              : 
     245           14 :          CALL generate_virtual_mo(qs_env, mo_set, evals_virt, mo_virt, nvirt, ispin=my_spin)
     246              :       ELSE
     247           46 :          NULLIFY (evals_virt)
     248           46 :          nvirt = 0
     249              :       END IF
     250              : 
     251           60 :       calc_matsh = .TRUE.
     252           60 :       IF (PRESENT(external_matrix_shalf)) calc_matsh = .FALSE.
     253              : 
     254              :       ! Create S^1/2 : from sparse to full matrix, if no external available
     255          116 :       IF (calc_matsh) THEN
     256           58 :          NULLIFY (matrix_shalf)
     257              :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     258           58 :                                   nrow_global=nrow_global, ncol_global=nrow_global)
     259           58 :          ALLOCATE (matrix_shalf)
     260           58 :          CALL cp_fm_create(matrix_shalf, fm_struct_tmp, name="matrix_shalf")
     261           58 :          CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_work")
     262           58 :          CALL cp_fm_struct_release(fm_struct_tmp)
     263           58 :          CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, matrix_shalf)
     264           58 :          CALL cp_fm_power(matrix_shalf, matrix_work, 0.5_dp, EPSILON(0.0_dp), n_dependent)
     265           58 :          CALL cp_fm_release(matrix_work)
     266              :       ELSE
     267              :          matrix_shalf => external_matrix_shalf
     268              :       END IF
     269              : 
     270              :       ! Multiply S^(1/2) time the mOS coefficients to get orthonormalized MOS
     271              :       CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     272           60 :                                nrow_global=nrow_global, ncol_global=ncol_global)
     273           60 :       CALL cp_fm_create(matrix_shalfc, fm_struct_tmp, name="matrix_shalfc")
     274              :       CALL parallel_gemm("N", "N", nrow_global, ncol_global, nrow_global, &
     275           60 :                          1.0_dp, matrix_shalf, mo_coeff, 0.0_dp, matrix_shalfc)
     276           60 :       CALL cp_fm_struct_release(fm_struct_tmp)
     277              : 
     278           60 :       IF (do_virt) THEN
     279           14 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T14,I10,T27,A))') &
     280            7 :             " Compute ", nvirt, " additional unoccupied KS orbitals"
     281              :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     282           14 :                                   nrow_global=nrow_global, ncol_global=nvirt)
     283           14 :          CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_shalfc")
     284              :          CALL parallel_gemm("N", "N", nrow_global, nvirt, nrow_global, &
     285           14 :                             1.0_dp, matrix_shalf, mo_virt, 0.0_dp, matrix_work)
     286           14 :          CALL cp_fm_struct_release(fm_struct_tmp)
     287              :       END IF
     288              : 
     289           60 :       IF (calc_matsh) THEN
     290           58 :          CALL cp_fm_release(matrix_shalf)
     291           58 :          DEALLOCATE (matrix_shalf)
     292              :       END IF
     293              :       ! Array to store the PDOS per kind and angular momentum
     294           60 :       do_ldos = .FALSE.
     295           60 :       ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%LDOS")
     296              : 
     297           60 :       CALL section_vals_get(ldos_section, n_repetition=nldos)
     298           60 :       IF (nldos > 0) THEN
     299            8 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
     300            4 :             " Prepare the list of atoms for LDOS.   Number of lists: ", nldos
     301            8 :          do_ldos = .TRUE.
     302           44 :          ALLOCATE (ldos_p(nldos))
     303           24 :          ALLOCATE (ldos_index(nldos))
     304           28 :          DO ildos = 1, nldos
     305           20 :             WRITE (ldos_index(ildos), '(I0)') ildos
     306           20 :             ALLOCATE (ldos_p(ildos)%ldos)
     307              :             NULLIFY (ldos_p(ildos)%ldos%pdos_array)
     308              :             NULLIFY (ldos_p(ildos)%ldos%list_index)
     309              : 
     310           20 :             CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep)
     311           20 :             IF (n_rep > 0) THEN
     312           20 :                ldos_p(ildos)%ldos%nlist = 0
     313           40 :                DO ir = 1, n_rep
     314           20 :                   NULLIFY (list)
     315              :                   CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, &
     316           20 :                                             i_vals=list)
     317           40 :                   IF (ASSOCIATED(list)) THEN
     318           20 :                      CALL reallocate(ldos_p(ildos)%ldos%list_index, 1, ldos_p(ildos)%ldos%nlist + SIZE(list))
     319           76 :                      DO i = 1, SIZE(list)
     320           76 :                         ldos_p(ildos)%ldos%list_index(i + ldos_p(ildos)%ldos%nlist) = list(i)
     321              :                      END DO
     322           20 :                      ldos_p(ildos)%ldos%nlist = ldos_p(ildos)%ldos%nlist + SIZE(list)
     323              :                   END IF
     324              :                END DO
     325              :             ELSE
     326              :                ! stop, LDOS without list of atoms is not implemented
     327              :             END IF
     328              : 
     329           20 :             IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') &
     330           10 :                " List ", ildos, " contains ", ldos_p(ildos)%ldos%nlist, " atoms"
     331              :             CALL section_vals_val_get(ldos_section, "COMPONENTS", i_rep_section=ildos, &
     332           20 :                                       l_val=ldos_p(ildos)%ldos%separate_components)
     333           20 :             IF (ldos_p(ildos)%ldos%separate_components) THEN
     334           16 :                ALLOCATE (ldos_p(ildos)%ldos%pdos_array(nsoset(maxlgto), nmo + nvirt))
     335              :             ELSE
     336           64 :                ALLOCATE (ldos_p(ildos)%ldos%pdos_array(0:maxlgto, nmo + nvirt))
     337              :             END IF
     338          716 :             ldos_p(ildos)%ldos%pdos_array = 0.0_dp
     339           48 :             ldos_p(ildos)%ldos%maxl = -1
     340              : 
     341              :          END DO
     342              :       END IF
     343              : 
     344           60 :       do_r_ldos = .FALSE.
     345           60 :       ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%R_LDOS")
     346           60 :       CALL section_vals_get(ldos_section, n_repetition=n_r_ldos)
     347           60 :       IF (n_r_ldos > 0) THEN
     348            0 :          do_r_ldos = .TRUE.
     349            0 :          IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') &
     350            0 :             " Prepare the list of points for R_LDOS.   Number of lists: ", n_r_ldos
     351            0 :          ALLOCATE (r_ldos_p(n_r_ldos))
     352            0 :          ALLOCATE (r_ldos_index(n_r_ldos))
     353              :          CALL get_qs_env(qs_env=qs_env, &
     354              :                          cell=cell, &
     355              :                          dft_control=dft_control, &
     356            0 :                          pw_env=pw_env)
     357              :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     358            0 :                          pw_pools=pw_pools)
     359              : 
     360            0 :          CALL auxbas_pw_pool%create_pw(wf_r)
     361            0 :          CALL auxbas_pw_pool%create_pw(wf_g)
     362            0 :          ALLOCATE (read_r(4, n_r_ldos))
     363            0 :          DO ildos = 1, n_r_ldos
     364            0 :             WRITE (r_ldos_index(ildos), '(I0)') ildos
     365            0 :             ALLOCATE (r_ldos_p(ildos)%ldos)
     366              :             NULLIFY (r_ldos_p(ildos)%ldos%pdos_array)
     367              :             NULLIFY (r_ldos_p(ildos)%ldos%list_index)
     368              : 
     369            0 :             CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep)
     370            0 :             IF (n_rep > 0) THEN
     371            0 :                r_ldos_p(ildos)%ldos%nlist = 0
     372            0 :                DO ir = 1, n_rep
     373            0 :                   NULLIFY (list)
     374              :                   CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, &
     375            0 :                                             i_vals=list)
     376            0 :                   IF (ASSOCIATED(list)) THEN
     377            0 :                      CALL reallocate(r_ldos_p(ildos)%ldos%list_index, 1, r_ldos_p(ildos)%ldos%nlist + SIZE(list))
     378            0 :                      DO i = 1, SIZE(list)
     379            0 :                         r_ldos_p(ildos)%ldos%list_index(i + r_ldos_p(ildos)%ldos%nlist) = list(i)
     380              :                      END DO
     381            0 :                      r_ldos_p(ildos)%ldos%nlist = r_ldos_p(ildos)%ldos%nlist + SIZE(list)
     382              :                   END IF
     383              :                END DO
     384              :             ELSE
     385              :                ! stop, LDOS without list of atoms is not implemented
     386              :             END IF
     387              : 
     388            0 :             ALLOCATE (r_ldos_p(ildos)%ldos%pdos_array(nmo + nvirt))
     389            0 :             r_ldos_p(ildos)%ldos%pdos_array = 0.0_dp
     390            0 :             read_r(1:3, ildos) = .FALSE.
     391            0 :             CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, explicit=read_r(1, ildos))
     392            0 :             IF (read_r(1, ildos)) THEN
     393              :                CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, r_vals= &
     394            0 :                                          r_ldos_p(ildos)%ldos%x_range)
     395              :             ELSE
     396            0 :                ALLOCATE (r_ldos_p(ildos)%ldos%x_range(2))
     397            0 :                r_ldos_p(ildos)%ldos%x_range = 0.0_dp
     398              :             END IF
     399            0 :             CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, explicit=read_r(2, ildos))
     400            0 :             IF (read_r(2, ildos)) THEN
     401              :                CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, r_vals= &
     402            0 :                                          r_ldos_p(ildos)%ldos%y_range)
     403              :             ELSE
     404            0 :                ALLOCATE (r_ldos_p(ildos)%ldos%y_range(2))
     405            0 :                r_ldos_p(ildos)%ldos%y_range = 0.0_dp
     406              :             END IF
     407            0 :             CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, explicit=read_r(3, ildos))
     408            0 :             IF (read_r(3, ildos)) THEN
     409              :                CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, r_vals= &
     410            0 :                                          r_ldos_p(ildos)%ldos%z_range)
     411              :             ELSE
     412            0 :                ALLOCATE (r_ldos_p(ildos)%ldos%z_range(2))
     413            0 :                r_ldos_p(ildos)%ldos%z_range = 0.0_dp
     414              :             END IF
     415              : 
     416            0 :             CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, explicit=read_r(4, ildos))
     417            0 :             IF (read_r(4, ildos)) THEN
     418              :                CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, &
     419            0 :                                          r_vals=r_ldos_p(ildos)%ldos%eval_range)
     420              :             ELSE
     421            0 :                ALLOCATE (r_ldos_p(ildos)%ldos%eval_range(2))
     422            0 :                r_ldos_p(ildos)%ldos%eval_range(1) = -HUGE(0.0_dp)
     423            0 :                r_ldos_p(ildos)%ldos%eval_range(2) = +HUGE(0.0_dp)
     424              :             END IF
     425              : 
     426            0 :             bo => wf_r%pw_grid%bounds_local
     427            0 :             dh = wf_r%pw_grid%dh
     428            0 :             dvol = wf_r%pw_grid%dvol
     429            0 :             np_tot = wf_r%pw_grid%npts(1)*wf_r%pw_grid%npts(2)*wf_r%pw_grid%npts(3)
     430            0 :             ALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local(3, np_tot))
     431              : 
     432            0 :             r_ldos_p(ildos)%ldos%npoints = 0
     433            0 :             DO jz = bo(1, 3), bo(2, 3)
     434            0 :             DO jy = bo(1, 2), bo(2, 2)
     435            0 :             DO jx = bo(1, 1), bo(2, 1)
     436              :                !compute the position of the grid point
     437            0 :                i = jx - wf_r%pw_grid%bounds(1, 1)
     438            0 :                j = jy - wf_r%pw_grid%bounds(1, 2)
     439            0 :                k = jz - wf_r%pw_grid%bounds(1, 3)
     440            0 :                r(3) = k*dh(3, 3) + j*dh(3, 2) + i*dh(3, 1)
     441            0 :                r(2) = k*dh(2, 3) + j*dh(2, 2) + i*dh(2, 1)
     442            0 :                r(1) = k*dh(1, 3) + j*dh(1, 2) + i*dh(1, 1)
     443              : 
     444            0 :                DO il = 1, r_ldos_p(ildos)%ldos%nlist
     445            0 :                   iatom = r_ldos_p(ildos)%ldos%list_index(il)
     446            0 :                   ratom = particle_set(iatom)%r
     447            0 :                   r_vec = pbc(ratom, r, cell)
     448            0 :                   IF (cell%orthorhombic) THEN
     449            0 :                      IF (cell%perd(1) == 0) r_vec(1) = MODULO(r_vec(1), cell%hmat(1, 1))
     450            0 :                      IF (cell%perd(2) == 0) r_vec(2) = MODULO(r_vec(2), cell%hmat(2, 2))
     451            0 :                      IF (cell%perd(3) == 0) r_vec(3) = MODULO(r_vec(3), cell%hmat(3, 3))
     452              :                   END IF
     453              : 
     454            0 :                   in_x = 0
     455            0 :                   in_y = 0
     456            0 :                   in_z = 0
     457            0 :                   IF (r_ldos_p(ildos)%ldos%x_range(1) /= 0.0_dp) THEN
     458            0 :                      IF (r_vec(1) > r_ldos_p(ildos)%ldos%x_range(1) .AND. &
     459              :                          r_vec(1) < r_ldos_p(ildos)%ldos%x_range(2)) THEN
     460            0 :                         in_x = 1
     461              :                      END IF
     462              :                   ELSE
     463              :                      in_x = 1
     464              :                   END IF
     465            0 :                   IF (r_ldos_p(ildos)%ldos%y_range(1) /= 0.0_dp) THEN
     466            0 :                      IF (r_vec(2) > r_ldos_p(ildos)%ldos%y_range(1) .AND. &
     467              :                          r_vec(2) < r_ldos_p(ildos)%ldos%y_range(2)) THEN
     468            0 :                         in_y = 1
     469              :                      END IF
     470              :                   ELSE
     471              :                      in_y = 1
     472              :                   END IF
     473            0 :                   IF (r_ldos_p(ildos)%ldos%z_range(1) /= 0.0_dp) THEN
     474            0 :                      IF (r_vec(3) > r_ldos_p(ildos)%ldos%z_range(1) .AND. &
     475              :                          r_vec(3) < r_ldos_p(ildos)%ldos%z_range(2)) THEN
     476            0 :                         in_z = 1
     477              :                      END IF
     478              :                   ELSE
     479              :                      in_z = 1
     480              :                   END IF
     481            0 :                   IF (in_x*in_y*in_z > 0) THEN
     482            0 :                      r_ldos_p(ildos)%ldos%npoints = r_ldos_p(ildos)%ldos%npoints + 1
     483            0 :                      r_ldos_p(ildos)%ldos%index_grid_local(1, r_ldos_p(ildos)%ldos%npoints) = jx
     484            0 :                      r_ldos_p(ildos)%ldos%index_grid_local(2, r_ldos_p(ildos)%ldos%npoints) = jy
     485            0 :                      r_ldos_p(ildos)%ldos%index_grid_local(3, r_ldos_p(ildos)%ldos%npoints) = jz
     486            0 :                      EXIT
     487              :                   END IF
     488              :                END DO
     489              :             END DO
     490              :             END DO
     491              :             END DO
     492            0 :             CALL reallocate(r_ldos_p(ildos)%ldos%index_grid_local, 1, 3, 1, r_ldos_p(ildos)%ldos%npoints)
     493            0 :             npoints = r_ldos_p(ildos)%ldos%npoints
     494            0 :             CALL para_env%sum(npoints)
     495            0 :             IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') &
     496            0 :                " List ", ildos, " contains ", npoints, " grid points"
     497              :          END DO
     498              :       END IF
     499              : 
     500           60 :       CALL section_vals_val_get(dft_section, "PRINT%PDOS%COMPONENTS", l_val=separate_components)
     501           60 :       IF (separate_components) THEN
     502           80 :          ALLOCATE (pdos_array(nsoset(maxlgto), nkind, nmo + nvirt))
     503              :       ELSE
     504          220 :          ALLOCATE (pdos_array(0:maxlgto, nkind, nmo + nvirt))
     505              :       END IF
     506           60 :       IF (do_virt) THEN
     507           42 :          ALLOCATE (eigenvalues(nmo + nvirt))
     508          122 :          eigenvalues(1:nmo) = mo_set%eigenvalues(1:nmo)
     509          346 :          eigenvalues(nmo + 1:nmo + nvirt) = evals_virt(1:nvirt)
     510           28 :          ALLOCATE (occupation_numbers(nmo + nvirt))
     511          288 :          occupation_numbers(:) = 0.0_dp
     512          122 :          occupation_numbers(1:nmo) = mo_set%occupation_numbers(1:nmo)
     513              :       ELSE
     514           46 :          eigenvalues => mo_set%eigenvalues
     515           46 :          occupation_numbers => mo_set%occupation_numbers
     516              :       END IF
     517              : 
     518         8036 :       pdos_array = 0.0_dp
     519           60 :       nao = mo_set%nao
     520          180 :       ALLOCATE (vecbuffer(1, nao))
     521         2432 :       vecbuffer = 0.0_dp
     522          180 :       ALLOCATE (firstrow(natom))
     523          224 :       firstrow = 0
     524              : 
     525              :       !Adjust energy range for r_ldos
     526           60 :       DO ildos = 1, n_r_ldos
     527            0 :          IF (eigenvalues(1) > r_ldos_p(ildos)%ldos%eval_range(1)) &
     528            0 :             r_ldos_p(ildos)%ldos%eval_range(1) = eigenvalues(1)
     529            0 :          IF (eigenvalues(nmo + nvirt) < r_ldos_p(ildos)%ldos%eval_range(2)) &
     530           60 :             r_ldos_p(ildos)%ldos%eval_range(2) = eigenvalues(nmo + nvirt)
     531              :       END DO
     532              : 
     533           60 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T15,A))') &
     534           30 :          "---- PDOS: start iteration on the KS states --- "
     535              : 
     536          548 :       DO imo = 1, nmo + nvirt
     537              : 
     538          488 :          IF (output_unit > 0 .AND. MOD(imo, out_each) == 0) WRITE (UNIT=output_unit, FMT='((T20,A,I10))') &
     539            0 :             " KS state index : ", imo
     540              :          ! Extract the eigenvector from the distributed full matrix
     541          488 :          IF (imo > nmo) THEN
     542              :             CALL cp_fm_get_submatrix(matrix_work, vecbuffer, 1, imo - nmo, &
     543          166 :                                      nao, 1, transpose=.TRUE.)
     544              :          ELSE
     545              :             CALL cp_fm_get_submatrix(matrix_shalfc, vecbuffer, 1, imo, &
     546          322 :                                      nao, 1, transpose=.TRUE.)
     547              :          END IF
     548              : 
     549              :          ! Calculate the pdos for all the kinds
     550          488 :          irow = 1
     551         1700 :          DO iatom = 1, natom
     552         1212 :             firstrow(iatom) = irow
     553         1212 :             NULLIFY (orb_basis_set)
     554         1212 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     555         1212 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     556              :             CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     557              :                                    nset=nset, &
     558              :                                    nshell=nshell, &
     559         1212 :                                    l=l, maxl=maxl)
     560         2912 :             IF (separate_components) THEN
     561          716 :                isgf = 1
     562         3256 :                DO iset = 1, nset
     563         6946 :                   DO ishell = 1, nshell(iset)
     564         3690 :                      lshell = l(ishell, iset)
     565        14544 :                      DO iso = 1, nso(lshell)
     566         8314 :                         lcomponent = nsoset(lshell - 1) + iso
     567              :                         pdos_array(lcomponent, ikind, imo) = &
     568              :                            pdos_array(lcomponent, ikind, imo) + &
     569         8314 :                            vecbuffer(1, irow)*vecbuffer(1, irow)
     570        12004 :                         irow = irow + 1
     571              :                      END DO ! iso
     572              :                   END DO ! ishell
     573              :                END DO ! iset
     574              :             ELSE
     575          496 :                isgf = 1
     576         1392 :                DO iset = 1, nset
     577         2928 :                   DO ishell = 1, nshell(iset)
     578         1536 :                      lshell = l(ishell, iset)
     579         5536 :                      DO iso = 1, nso(lshell)
     580              :                         pdos_array(lshell, ikind, imo) = &
     581              :                            pdos_array(lshell, ikind, imo) + &
     582         3104 :                            vecbuffer(1, irow)*vecbuffer(1, irow)
     583         4640 :                         irow = irow + 1
     584              :                      END DO ! iso
     585              :                   END DO ! ishell
     586              :                END DO ! iset
     587              :             END IF
     588              :          END DO ! iatom
     589              : 
     590              :          ! Calculate the pdos for all the lists
     591          608 :          DO ildos = 1, nldos
     592          932 :             DO il = 1, ldos_p(ildos)%ldos%nlist
     593          324 :                iatom = ldos_p(ildos)%ldos%list_index(il)
     594              : 
     595          324 :                irow = firstrow(iatom)
     596          324 :                NULLIFY (orb_basis_set)
     597          324 :                CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     598          324 :                CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     599              : 
     600              :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     601              :                                       nset=nset, &
     602              :                                       nshell=nshell, &
     603          324 :                                       l=l, maxl=maxl)
     604          324 :                ldos_p(ildos)%ldos%maxl = MAX(ldos_p(ildos)%ldos%maxl, maxl)
     605          768 :                IF (ldos_p(ildos)%ldos%separate_components) THEN
     606          108 :                   isgf = 1
     607          324 :                   DO iset = 1, nset
     608          720 :                      DO ishell = 1, nshell(iset)
     609          396 :                         lshell = l(ishell, iset)
     610         1440 :                         DO iso = 1, nso(lshell)
     611          828 :                            lcomponent = nsoset(lshell - 1) + iso
     612              :                            ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) = &
     613              :                               ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) + &
     614          828 :                               vecbuffer(1, irow)*vecbuffer(1, irow)
     615         1224 :                            irow = irow + 1
     616              :                         END DO ! iso
     617              :                      END DO ! ishell
     618              :                   END DO ! iset
     619              :                ELSE
     620          216 :                   isgf = 1
     621          648 :                   DO iset = 1, nset
     622         1428 :                      DO ishell = 1, nshell(iset)
     623          780 :                         lshell = l(ishell, iset)
     624         2820 :                         DO iso = 1, nso(lshell)
     625              :                            ldos_p(ildos)%ldos%pdos_array(lshell, imo) = &
     626              :                               ldos_p(ildos)%ldos%pdos_array(lshell, imo) + &
     627         1608 :                               vecbuffer(1, irow)*vecbuffer(1, irow)
     628         2388 :                            irow = irow + 1
     629              :                         END DO ! iso
     630              :                      END DO ! ishell
     631              :                   END DO ! iset
     632              :                END IF
     633              :             END DO !il
     634              :          END DO !ildos
     635              : 
     636              :          ! Calculate the DOS projected in a given volume in real space
     637          548 :          DO ildos = 1, n_r_ldos
     638            0 :             IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. &
     639          488 :                 r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN
     640              : 
     641            0 :                IF (imo > nmo) THEN
     642              :                   CALL calculate_wavefunction(mo_virt, imo - nmo, &
     643              :                                               wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
     644            0 :                                               pw_env)
     645              :                ELSE
     646              :                   CALL calculate_wavefunction(mo_coeff, imo, &
     647              :                                               wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
     648            0 :                                               pw_env)
     649              :                END IF
     650            0 :                r_ldos_p(ildos)%ldos%pdos_array(imo) = 0.0_dp
     651            0 :                DO il = 1, r_ldos_p(ildos)%ldos%npoints
     652            0 :                   j = j + 1
     653            0 :                   jx = r_ldos_p(ildos)%ldos%index_grid_local(1, il)
     654            0 :                   jy = r_ldos_p(ildos)%ldos%index_grid_local(2, il)
     655            0 :                   jz = r_ldos_p(ildos)%ldos%index_grid_local(3, il)
     656              :                   r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo) + &
     657            0 :                                                          wf_r%array(jx, jy, jz)*wf_r%array(jx, jy, jz)
     658              :                END DO
     659            0 :                r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo)*dvol
     660              :             END IF
     661              :          END DO
     662              :       END DO ! imo
     663              : 
     664           60 :       CALL cp_fm_release(matrix_shalfc)
     665           60 :       DEALLOCATE (vecbuffer)
     666              : 
     667           60 :       CALL section_vals_val_get(dft_section, "PRINT%PDOS%APPEND", l_val=append)
     668           60 :       IF (append .AND. iterstep > 1) THEN
     669            6 :          my_pos = "APPEND"
     670              :       ELSE
     671           54 :          my_pos = "REWIND"
     672              :       END IF
     673           60 :       my_act = "WRITE"
     674          176 :       DO ikind = 1, nkind
     675              : 
     676          116 :          NULLIFY (orb_basis_set)
     677          116 :          CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name)
     678          116 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     679          116 :          CALL get_gto_basis_set(gto_basis_set=orb_basis_set, maxl=maxl)
     680              : 
     681              :          ! basis none has no associated maxl, and no pdos
     682          116 :          IF (maxl < 0) CYCLE
     683              : 
     684          116 :          IF (PRESENT(ispin)) THEN
     685           20 :             IF (PRESENT(xas_mittle)) THEN
     686           20 :                my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
     687              :             ELSE
     688            0 :                my_mittle = TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind)))
     689              :             END IF
     690           20 :             my_spin = ispin
     691              :          ELSE
     692           96 :             my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind)))
     693           96 :             my_spin = 1
     694              :          END IF
     695              : 
     696              :          iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
     697              :                                    extension=".pdos", file_position=my_pos, file_action=my_act, &
     698          116 :                                    file_form="FORMATTED", middle_name=TRIM(my_mittle))
     699          116 :          IF (iw > 0) THEN
     700              : 
     701           58 :             fmtstr1 = "(I8,2X,2F16.6,  (2X,F16.8))"
     702           58 :             fmtstr2 = "(A42,  (10X,A8))"
     703           58 :             IF (separate_components) THEN
     704           16 :                WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(maxl)
     705           16 :                WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(maxl)
     706              :             ELSE
     707           42 :                WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") maxl + 1
     708           42 :                WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") maxl + 1
     709              :             END IF
     710              : 
     711              :             WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,A)") &
     712           58 :                "# Projected DOS for atomic kind "//TRIM(kind_name)//" at iteration step i = ", &
     713          116 :                iterstep, ", E(Fermi) = ", e_fermi, " a.u."
     714           58 :             IF (separate_components) THEN
     715           64 :                ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl))
     716          496 :                tmp_str = ""
     717           60 :                DO j = 0, maxl
     718          184 :                   DO i = -j, j
     719          168 :                      tmp_str(0, j, i) = sgf_symbol(0, j, i)
     720              :                   END DO
     721              :                END DO
     722              : 
     723              :                WRITE (UNIT=iw, FMT=fmtstr2) &
     724           16 :                   "#     MO Eigenvalue [a.u.]      Occupation", &
     725          200 :                   ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, maxl)
     726          328 :                DO imo = 1, nmo + nvirt
     727          312 :                   WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
     728          640 :                      (pdos_array(lshell, ikind, imo), lshell=1, nsoset(maxl))
     729              :                END DO
     730           16 :                DEALLOCATE (tmp_str)
     731              :             ELSE
     732              :                WRITE (UNIT=iw, FMT=fmtstr2) &
     733           42 :                   "#     MO Eigenvalue [a.u.]      Occupation", &
     734          178 :                   (TRIM(l_sym(il)), il=0, maxl)
     735          210 :                DO imo = 1, nmo + nvirt
     736          168 :                   WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
     737          378 :                      (pdos_array(lshell, ikind, imo), lshell=0, maxl)
     738              :                END DO
     739              :             END IF
     740              :          END IF
     741              :          CALL cp_print_key_finished_output(iw, logger, dft_section, &
     742          292 :                                            "PRINT%PDOS")
     743              : 
     744              :       END DO ! ikind
     745              : 
     746              :       ! write the pdos for the lists, each ona different file,
     747              :       ! the filenames are indexed with the list number
     748           80 :       DO ildos = 1, nldos
     749              :          ! basis none has no associated maxl, and no pdos
     750           80 :          IF (ldos_p(ildos)%ldos%maxl > 0) THEN
     751              : 
     752           20 :             IF (PRESENT(ispin)) THEN
     753            0 :                IF (PRESENT(xas_mittle)) THEN
     754            0 :                   my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos))
     755              :                ELSE
     756            0 :                   my_mittle = TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos))
     757              :                END IF
     758            0 :                my_spin = ispin
     759              :             ELSE
     760           20 :                my_mittle = "list"//TRIM(ldos_index(ildos))
     761           20 :                my_spin = 1
     762              :             END IF
     763              : 
     764              :             iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
     765              :                                       extension=".pdos", file_position=my_pos, file_action=my_act, &
     766           20 :                                       file_form="FORMATTED", middle_name=TRIM(my_mittle))
     767           20 :             IF (iw > 0) THEN
     768              : 
     769           10 :                fmtstr1 = "(I8,2X,2F16.6,  (2X,F16.8))"
     770           10 :                fmtstr2 = "(A42,  (10X,A8))"
     771           10 :                IF (ldos_p(ildos)%ldos%separate_components) THEN
     772            2 :                   WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl)
     773            2 :                   WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl)
     774              :                ELSE
     775            8 :                   WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 1
     776            8 :                   WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 1
     777              :                END IF
     778              : 
     779              :                WRITE (UNIT=iw, FMT="(A,I0,A,I0,A,I0,A,F12.6,A)") &
     780           10 :                   "# Projected DOS for list ", ildos, " of ", ldos_p(ildos)%ldos%nlist, &
     781           10 :                   " atoms, at iteration step i = ", iterstep, &
     782           20 :                   ", E(Fermi) = ", e_fermi, " a.u."
     783           10 :                IF (ldos_p(ildos)%ldos%separate_components) THEN
     784            8 :                   ALLOCATE (tmp_str(0:0, 0:ldos_p(ildos)%ldos%maxl, -ldos_p(ildos)%ldos%maxl:ldos_p(ildos)%ldos%maxl))
     785           72 :                   tmp_str = ""
     786            8 :                   DO j = 0, ldos_p(ildos)%ldos%maxl
     787           26 :                      DO i = -j, j
     788           24 :                         tmp_str(0, j, i) = sgf_symbol(0, j, i)
     789              :                      END DO
     790              :                   END DO
     791              : 
     792              :                   WRITE (UNIT=iw, FMT=fmtstr2) &
     793            2 :                      "#     MO Eigenvalue [a.u.]      Occupation", &
     794           28 :                      ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, ldos_p(ildos)%ldos%maxl)
     795           20 :                   DO imo = 1, nmo + nvirt
     796           18 :                      WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
     797          200 :                         (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=1, nsoset(ldos_p(ildos)%ldos%maxl))
     798              :                   END DO
     799            2 :                   DEALLOCATE (tmp_str)
     800              :                ELSE
     801              :                   WRITE (UNIT=iw, FMT=fmtstr2) &
     802            8 :                      "#     MO Eigenvalue [a.u.]      Occupation", &
     803           39 :                      (TRIM(l_sym(il)), il=0, ldos_p(ildos)%ldos%maxl)
     804           50 :                   DO imo = 1, nmo + nvirt
     805           42 :                      WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), &
     806          209 :                         (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=0, ldos_p(ildos)%ldos%maxl)
     807              :                   END DO
     808              :                END IF
     809              :             END IF
     810              :             CALL cp_print_key_finished_output(iw, logger, dft_section, &
     811           20 :                                               "PRINT%PDOS")
     812              :          END IF ! maxl>0
     813              :       END DO ! ildos
     814              : 
     815              :       ! write the pdos for the lists, each ona different file,
     816              :       ! the filenames are indexed with the list number
     817           60 :       DO ildos = 1, n_r_ldos
     818              : 
     819            0 :          npoints = r_ldos_p(ildos)%ldos%npoints
     820            0 :          CALL para_env%sum(npoints)
     821            0 :          CALL para_env%sum(np_tot)
     822            0 :          CALL para_env%sum(r_ldos_p(ildos)%ldos%pdos_array)
     823            0 :          IF (PRESENT(ispin)) THEN
     824            0 :             IF (PRESENT(xas_mittle)) THEN
     825            0 :                my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos))
     826              :             ELSE
     827            0 :                my_mittle = TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos))
     828              :             END IF
     829            0 :             my_spin = ispin
     830              :          ELSE
     831            0 :             my_mittle = "r_list"//TRIM(r_ldos_index(ildos))
     832            0 :             my_spin = 1
     833              :          END IF
     834              : 
     835              :          iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", &
     836              :                                    extension=".pdos", file_position=my_pos, file_action=my_act, &
     837            0 :                                    file_form="FORMATTED", middle_name=TRIM(my_mittle))
     838            0 :          IF (iw > 0) THEN
     839            0 :             fmtstr1 = "(I8,2X,2F16.6,  (2X,F16.8))"
     840            0 :             fmtstr2 = "(A42,  (10X,A8))"
     841              : 
     842              :             WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,F12.6,A,F12.6,A)") &
     843            0 :                "# Projected DOS in real space, using ", npoints, &
     844            0 :                " points of the grid, and eval in the range", r_ldos_p(ildos)%ldos%eval_range(1:2), &
     845            0 :                " Hartree, E(Fermi) = ", e_fermi, " a.u."
     846              :             WRITE (UNIT=iw, FMT="(A)") &
     847            0 :                "#     MO Eigenvalue [a.u.]      Occupation      LDOS"
     848            0 :             DO imo = 1, nmo + nvirt
     849            0 :                IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. &
     850            0 :                    r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN
     851            0 :                   WRITE (UNIT=iw, FMT="(I8,2X,2F16.6,E20.10,E20.10)") imo, eigenvalues(imo), occupation_numbers(imo), &
     852            0 :                      r_ldos_p(ildos)%ldos%pdos_array(imo), r_ldos_p(ildos)%ldos%pdos_array(imo)*np_tot
     853              :                END IF
     854              :             END DO
     855              : 
     856              :          END IF
     857              :          CALL cp_print_key_finished_output(iw, logger, dft_section, &
     858           60 :                                            "PRINT%PDOS")
     859              :       END DO
     860              : 
     861              :       ! deallocate local variables
     862           60 :       DEALLOCATE (pdos_array)
     863           60 :       DEALLOCATE (firstrow)
     864           60 :       IF (do_ldos) THEN
     865           28 :          DO ildos = 1, nldos
     866           20 :             DEALLOCATE (ldos_p(ildos)%ldos%pdos_array)
     867           20 :             DEALLOCATE (ldos_p(ildos)%ldos%list_index)
     868           28 :             DEALLOCATE (ldos_p(ildos)%ldos)
     869              :          END DO
     870            8 :          DEALLOCATE (ldos_p)
     871            8 :          DEALLOCATE (ldos_index)
     872              :       END IF
     873           60 :       IF (do_r_ldos) THEN
     874            0 :          DO ildos = 1, n_r_ldos
     875            0 :             DEALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local)
     876            0 :             DEALLOCATE (r_ldos_p(ildos)%ldos%pdos_array)
     877            0 :             DEALLOCATE (r_ldos_p(ildos)%ldos%list_index)
     878            0 :             IF (.NOT. read_r(1, ildos)) THEN
     879            0 :                DEALLOCATE (r_ldos_p(ildos)%ldos%x_range)
     880              :             END IF
     881            0 :             IF (.NOT. read_r(2, ildos)) THEN
     882            0 :                DEALLOCATE (r_ldos_p(ildos)%ldos%y_range)
     883              :             END IF
     884            0 :             IF (.NOT. read_r(3, ildos)) THEN
     885            0 :                DEALLOCATE (r_ldos_p(ildos)%ldos%z_range)
     886              :             END IF
     887            0 :             IF (.NOT. read_r(4, ildos)) THEN
     888            0 :                DEALLOCATE (r_ldos_p(ildos)%ldos%eval_range)
     889              :             END IF
     890            0 :             DEALLOCATE (r_ldos_p(ildos)%ldos)
     891              :          END DO
     892            0 :          DEALLOCATE (read_r)
     893            0 :          DEALLOCATE (r_ldos_p)
     894            0 :          DEALLOCATE (r_ldos_index)
     895            0 :          CALL auxbas_pw_pool%give_back_pw(wf_r)
     896            0 :          CALL auxbas_pw_pool%give_back_pw(wf_g)
     897              :       END IF
     898           60 :       IF (do_virt) THEN
     899           14 :          DEALLOCATE (evals_virt)
     900           14 :          CALL cp_fm_release(mo_virt)
     901           14 :          CALL cp_fm_release(matrix_work)
     902           14 :          DEALLOCATE (eigenvalues)
     903           14 :          DEALLOCATE (occupation_numbers)
     904              :       END IF
     905              : 
     906           60 :       CALL timestop(handle)
     907              : 
     908          600 :    END SUBROUTINE calculate_projected_dos
     909              : 
     910              : ! **************************************************************************************************
     911              : !> \brief   Compute additional virtual states  starting from the available MOS
     912              : !> \param qs_env ...
     913              : !> \param mo_set ...
     914              : !> \param evals_virt ...
     915              : !> \param mo_virt ...
     916              : !> \param nvirt ...
     917              : !> \param ispin ...
     918              : !> \date    08.03.2008
     919              : !> \par History:
     920              : !>       -
     921              : !> \par Variables
     922              : !>       -
     923              : !>       -
     924              : !> \author  MI
     925              : !> \version 1.0
     926              : ! **************************************************************************************************
     927              : 
     928           14 :    SUBROUTINE generate_virtual_mo(qs_env, mo_set, evals_virt, mo_virt, &
     929              :                                   nvirt, ispin)
     930              : 
     931              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     932              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     933              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: evals_virt
     934              :       TYPE(cp_fm_type), INTENT(OUT)                      :: mo_virt
     935              :       INTEGER, INTENT(IN)                                :: nvirt, ispin
     936              : 
     937              :       INTEGER                                            :: nmo, nrow_global
     938              :       TYPE(cp_blacs_env_type), POINTER                   :: context
     939              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     940              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     941           14 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix, s_matrix
     942              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     943              :       TYPE(preconditioner_type), POINTER                 :: local_preconditioner
     944              :       TYPE(scf_control_type), POINTER                    :: scf_control
     945              : 
     946           14 :       NULLIFY (evals_virt)
     947           42 :       ALLOCATE (evals_virt(nvirt))
     948              : 
     949              :       CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, &
     950           14 :                       scf_control=scf_control)
     951           14 :       CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, nmo=nmo)
     952              :       CALL cp_fm_get_info(mo_coeff, context=context, para_env=para_env, &
     953           14 :                           nrow_global=nrow_global)
     954              : 
     955              :       CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     956           14 :                                nrow_global=nrow_global, ncol_global=nvirt)
     957           14 :       CALL cp_fm_create(mo_virt, fm_struct_tmp, name="virtual")
     958           14 :       CALL cp_fm_struct_release(fm_struct_tmp)
     959           14 :       CALL cp_fm_init_random(mo_virt, nvirt)
     960              : 
     961           14 :       NULLIFY (local_preconditioner)
     962              : 
     963              :       CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
     964              :                           matrix_c_fm=mo_virt, matrix_orthogonal_space_fm=mo_coeff, &
     965              :                           eps_gradient=scf_control%eps_lumos, &
     966              :                           preconditioner=local_preconditioner, &
     967              :                           iter_max=scf_control%max_iter_lumos, &
     968           14 :                           size_ortho_space=nmo)
     969              : 
     970              :       CALL calculate_subspace_eigenvalues(mo_virt, ks_matrix(ispin)%matrix, &
     971           14 :                                           evals_virt)
     972              : 
     973           42 :    END SUBROUTINE generate_virtual_mo
     974              : 
     975            0 : END MODULE qs_pdos
     976              : 
        

Generated by: LCOV version 2.0-1