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

            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 Some utilities for the construction of
      10              : !>      the localization environment
      11              : !> \author MI (05-2005)
      12              : ! **************************************************************************************************
      13              : MODULE qs_loc_utils
      14              : 
      15              :    USE ai_moments,                      ONLY: contract_cossin,&
      16              :                                               cossin
      17              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      18              :                                               gto_basis_set_type
      19              :    USE block_p_types,                   ONLY: block_p_type
      20              :    USE cell_types,                      ONLY: cell_type,&
      21              :                                               pbc
      22              :    USE cp_array_utils,                  ONLY: cp_1d_r_p_type
      23              :    USE cp_control_types,                ONLY: dft_control_type
      24              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      25              :                                               dbcsr_get_block_p,&
      26              :                                               dbcsr_p_type,&
      27              :                                               dbcsr_set,&
      28              :                                               dbcsr_type
      29              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      30              :    USE cp_files,                        ONLY: close_file,&
      31              :                                               open_file
      32              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale
      33              :    USE cp_fm_diag,                      ONLY: choose_eigv_solver
      34              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      35              :                                               cp_fm_struct_release,&
      36              :                                               cp_fm_struct_type
      37              :    USE cp_fm_types,                     ONLY: &
      38              :         cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
      39              :         cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type, cp_fm_write_unformatted
      40              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      41              :                                               cp_logger_get_default_io_unit,&
      42              :                                               cp_logger_type,&
      43              :                                               cp_to_string
      44              :    USE cp_output_handling,              ONLY: cp_p_file,&
      45              :                                               cp_print_key_finished_output,&
      46              :                                               cp_print_key_generate_filename,&
      47              :                                               cp_print_key_should_output,&
      48              :                                               cp_print_key_unit_nr
      49              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      50              :    USE input_constants,                 ONLY: &
      51              :         do_loc_crazy, do_loc_direct, do_loc_gapo, do_loc_jacobi, do_loc_l1_norm_sd, do_loc_none, &
      52              :         do_loc_scdm, energy_loc_range, op_loc_berry, op_loc_boys, op_loc_pipek, state_loc_all, &
      53              :         state_loc_list, state_loc_mixed, state_loc_none, state_loc_range
      54              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      55              :                                               section_vals_type,&
      56              :                                               section_vals_val_get
      57              :    USE kinds,                           ONLY: default_path_length,&
      58              :                                               default_string_length,&
      59              :                                               dp
      60              :    USE mathconstants,                   ONLY: twopi
      61              :    USE memory_utilities,                ONLY: reallocate
      62              :    USE message_passing,                 ONLY: mp_para_env_type
      63              :    USE orbital_pointers,                ONLY: ncoset
      64              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      65              :    USE particle_types,                  ONLY: particle_type
      66              :    USE qs_environment_types,            ONLY: get_qs_env,&
      67              :                                               qs_environment_type
      68              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      69              :                                               get_qs_kind_set,&
      70              :                                               qs_kind_type
      71              :    USE qs_loc_types,                    ONLY: get_qs_loc_env,&
      72              :                                               localized_wfn_control_create,&
      73              :                                               localized_wfn_control_release,&
      74              :                                               localized_wfn_control_type,&
      75              :                                               qs_loc_env_type,&
      76              :                                               set_qs_loc_env
      77              :    USE qs_localization_methods,         ONLY: initialize_weights
      78              :    USE qs_mo_methods,                   ONLY: make_mo_eig
      79              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      80              :                                               mo_set_type
      81              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      82              :                                               neighbor_list_iterate,&
      83              :                                               neighbor_list_iterator_create,&
      84              :                                               neighbor_list_iterator_p_type,&
      85              :                                               neighbor_list_iterator_release,&
      86              :                                               neighbor_list_set_p_type
      87              :    USE qs_scf_types,                    ONLY: ot_method_nr
      88              :    USE scf_control_types,               ONLY: scf_control_type
      89              : #include "./base/base_uses.f90"
      90              : 
      91              :    IMPLICIT NONE
      92              : 
      93              :    PRIVATE
      94              : 
      95              : ! *** Global parameters ***
      96              : 
      97              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_utils'
      98              : 
      99              : ! *** Public ***
     100              :    PUBLIC :: qs_loc_env_init, loc_write_restart, &
     101              :              retain_history, qs_loc_init, compute_berry_operator, &
     102              :              set_loc_centers, set_loc_wfn_lists, qs_loc_control_init
     103              : 
     104              : CONTAINS
     105              : 
     106              : ! **************************************************************************************************
     107              : !> \brief copy old mos to new ones, allocating as necessary
     108              : !> \param mo_loc_history ...
     109              : !> \param mo_loc ...
     110              : ! **************************************************************************************************
     111           10 :    SUBROUTINE retain_history(mo_loc_history, mo_loc)
     112              : 
     113              :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mo_loc_history
     114              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_loc
     115              : 
     116              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'retain_history'
     117              : 
     118              :       INTEGER                                            :: handle, i, ncol_hist, ncol_loc
     119              : 
     120           10 :       CALL timeset(routineN, handle)
     121              : 
     122           10 :       IF (.NOT. ASSOCIATED(mo_loc_history)) THEN
     123            8 :          ALLOCATE (mo_loc_history(SIZE(mo_loc)))
     124            4 :          DO i = 1, SIZE(mo_loc_history)
     125            4 :             CALL cp_fm_create(mo_loc_history(i), mo_loc(i)%matrix_struct)
     126              :          END DO
     127              :       END IF
     128              : 
     129           20 :       DO i = 1, SIZE(mo_loc_history)
     130           10 :          CALL cp_fm_get_info(mo_loc_history(i), ncol_global=ncol_hist)
     131           10 :          CALL cp_fm_get_info(mo_loc(i), ncol_global=ncol_loc)
     132           10 :          CPASSERT(ncol_hist == ncol_loc)
     133           30 :          CALL cp_fm_to_fm(mo_loc(i), mo_loc_history(i))
     134              :       END DO
     135              : 
     136           10 :       CALL timestop(handle)
     137              : 
     138           10 :    END SUBROUTINE retain_history
     139              : 
     140              : ! **************************************************************************************************
     141              : !> \brief rotate the mo_new, so that the orbitals are as similar
     142              : !>        as possible to ones in mo_ref.
     143              : !> \param mo_new ...
     144              : !> \param mo_ref ...
     145              : !> \param matrix_S ...
     146              : ! **************************************************************************************************
     147            8 :    SUBROUTINE rotate_state_to_ref(mo_new, mo_ref, matrix_S)
     148              : 
     149              :       TYPE(cp_fm_type), INTENT(IN)                       :: mo_new, mo_ref
     150              :       TYPE(dbcsr_type), POINTER                          :: matrix_S
     151              : 
     152              :       CHARACTER(len=*), PARAMETER :: routineN = 'rotate_state_to_ref'
     153              : 
     154              :       INTEGER                                            :: handle, ncol, ncol_ref, nrow
     155              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues
     156              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     157              :       TYPE(cp_fm_type)                                   :: o1, o2, o3, o4, smo
     158              : 
     159            8 :       CALL timeset(routineN, handle)
     160              : 
     161            8 :       CALL cp_fm_get_info(mo_new, nrow_global=nrow, ncol_global=ncol)
     162            8 :       CALL cp_fm_get_info(mo_ref, ncol_global=ncol_ref)
     163            8 :       CPASSERT(ncol == ncol_ref)
     164              : 
     165            8 :       NULLIFY (fm_struct_tmp)
     166            8 :       CALL cp_fm_create(smo, mo_ref%matrix_struct)
     167              : 
     168              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, &
     169              :                                ncol_global=ncol, para_env=mo_new%matrix_struct%para_env, &
     170            8 :                                context=mo_new%matrix_struct%context)
     171            8 :       CALL cp_fm_create(o1, fm_struct_tmp)
     172            8 :       CALL cp_fm_create(o2, fm_struct_tmp)
     173            8 :       CALL cp_fm_create(o3, fm_struct_tmp)
     174            8 :       CALL cp_fm_create(o4, fm_struct_tmp)
     175            8 :       CALL cp_fm_struct_release(fm_struct_tmp)
     176              : 
     177              :       ! o1 = (mo_new)^T matrix_S mo_ref
     178            8 :       CALL cp_dbcsr_sm_fm_multiply(matrix_S, mo_ref, smo, ncol)
     179            8 :       CALL parallel_gemm('T', 'N', ncol, ncol, nrow, 1.0_dp, mo_new, smo, 0.0_dp, o1)
     180              : 
     181              :       ! o2 = (o1^T o1)
     182            8 :       CALL parallel_gemm('T', 'N', ncol, ncol, ncol, 1.0_dp, o1, o1, 0.0_dp, o2)
     183              : 
     184              :       ! o2 = (o1^T o1)^-1/2
     185           24 :       ALLOCATE (eigenvalues(ncol))
     186            8 :       CALL choose_eigv_solver(o2, o3, eigenvalues)
     187            8 :       CALL cp_fm_to_fm(o3, o4)
     188           72 :       eigenvalues(:) = 1.0_dp/SQRT(eigenvalues(:))
     189            8 :       CALL cp_fm_column_scale(o4, eigenvalues)
     190            8 :       CALL parallel_gemm('N', 'T', ncol, ncol, ncol, 1.0_dp, o3, o4, 0.0_dp, o2)
     191              : 
     192              :       ! o3 = o1 (o1^T o1)^-1/2
     193            8 :       CALL parallel_gemm('N', 'N', ncol, ncol, ncol, 1.0_dp, o1, o2, 0.0_dp, o3)
     194              : 
     195              :       ! mo_new o1 (o1^T o1)^-1/2
     196            8 :       CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, mo_new, o3, 0.0_dp, smo)
     197            8 :       CALL cp_fm_to_fm(smo, mo_new)
     198              : 
     199              :       ! XXXXXXX testing
     200              :       ! CALL parallel_gemm('N','T',ncol,ncol,ncol,1.0_dp,o3,o3,0.0_dp,o1)
     201              :       ! WRITE(*,*) o1%local_data
     202              :       ! CALL parallel_gemm('T','N',ncol,ncol,ncol,1.0_dp,o3,o3,0.0_dp,o1)
     203              :       ! WRITE(*,*) o1%local_data
     204              : 
     205            8 :       CALL cp_fm_release(o1)
     206            8 :       CALL cp_fm_release(o2)
     207            8 :       CALL cp_fm_release(o3)
     208            8 :       CALL cp_fm_release(o4)
     209            8 :       CALL cp_fm_release(smo)
     210              : 
     211            8 :       CALL timestop(handle)
     212              : 
     213           32 :    END SUBROUTINE rotate_state_to_ref
     214              : 
     215              : ! **************************************************************************************************
     216              : !> \brief allocates the data, and initializes the operators
     217              : !> \param qs_loc_env new environment for the localization calculations
     218              : !> \param localized_wfn_control variables and directives for the localization
     219              : !> \param qs_env the qs_env in which the qs_env lives
     220              : !> \param myspin ...
     221              : !> \param do_localize ...
     222              : !> \param loc_coeff ...
     223              : !> \param mo_loc_history ...
     224              : !> \par History
     225              : !>      04.2005 created [MI]
     226              : !> \author MI
     227              : !> \note
     228              : !>      similar to the old one, but not quite
     229              : ! **************************************************************************************************
     230          892 :    SUBROUTINE qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, &
     231          446 :                               loc_coeff, mo_loc_history)
     232              : 
     233              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
     234              :       TYPE(localized_wfn_control_type), POINTER          :: localized_wfn_control
     235              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     236              :       INTEGER, INTENT(IN), OPTIONAL                      :: myspin
     237              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_localize
     238              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
     239              :          OPTIONAL                                        :: loc_coeff
     240              :       TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER  :: mo_loc_history
     241              : 
     242              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_loc_env_init'
     243              : 
     244              :       INTEGER                                            :: dim_op, handle, i, iatom, imo, imoloc, &
     245              :                                                             ispin, j, l_spin, lb, nao, naosub, &
     246              :                                                             natoms, nmo, nmosub, nspins, s_spin, ub
     247              :       REAL(KIND=dp)                                      :: my_occ, occ_imo
     248          446 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: occupations
     249          446 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
     250              :       TYPE(cell_type), POINTER                           :: cell
     251              :       TYPE(cp_fm_struct_type), POINTER                   :: tmp_fm_struct
     252          446 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: moloc_coeff
     253              :       TYPE(cp_fm_type), POINTER                          :: mat_ptr, mo_coeff
     254          446 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     255              :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     256          446 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     257              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     258          446 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     259              : 
     260          446 :       CALL timeset(routineN, handle)
     261              : 
     262          446 :       NULLIFY (mos, matrix_s, moloc_coeff, particle_set, para_env, cell, local_molecules, occupations, mat_ptr)
     263          446 :       IF (PRESENT(do_localize)) qs_loc_env%do_localize = do_localize
     264          446 :       IF (qs_loc_env%do_localize) THEN
     265              :          CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, cell=cell, &
     266              :                          local_molecules=local_molecules, particle_set=particle_set, &
     267          446 :                          para_env=para_env, mos=mos)
     268          446 :          nspins = SIZE(mos, 1)
     269          446 :          s_spin = 1
     270          446 :          l_spin = nspins
     271          446 :          IF (PRESENT(myspin)) THEN
     272          158 :             s_spin = myspin
     273          158 :             l_spin = myspin
     274              :          END IF
     275         1896 :          ALLOCATE (moloc_coeff(s_spin:l_spin))
     276         1004 :          DO ispin = s_spin, l_spin
     277          558 :             NULLIFY (tmp_fm_struct, mo_coeff)
     278          558 :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
     279          558 :             nmosub = localized_wfn_control%nloc_states(ispin)
     280              :             CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
     281          558 :                                      ncol_global=nmosub, para_env=para_env, context=mo_coeff%matrix_struct%context)
     282          558 :             CALL cp_fm_create(moloc_coeff(ispin), tmp_fm_struct)
     283              : 
     284              :             CALL cp_fm_get_info(moloc_coeff(ispin), nrow_global=naosub, &
     285          558 :                                 ncol_global=nmosub)
     286          558 :             CPASSERT(nao == naosub)
     287          558 :             IF ((localized_wfn_control%do_homo) .OR. &
     288              :                 (localized_wfn_control%set_of_states == state_loc_mixed)) THEN
     289          546 :                CPASSERT(nmo >= nmosub)
     290              :             ELSE
     291           12 :                CPASSERT(nao - nmo >= nmosub)
     292              :             END IF
     293          558 :             CALL cp_fm_set_all(moloc_coeff(ispin), 0.0_dp)
     294         2120 :             CALL cp_fm_struct_release(tmp_fm_struct)
     295              :          END DO ! ispin
     296              :          ! Copy the submatrix
     297              : 
     298          446 :          IF (PRESENT(loc_coeff)) ALLOCATE (mat_ptr)
     299              : 
     300         1004 :          DO ispin = s_spin, l_spin
     301              :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
     302          558 :                             occupation_numbers=occupations, nao=nao, nmo=nmo)
     303          558 :             lb = localized_wfn_control%lu_bound_states(1, ispin)
     304          558 :             ub = localized_wfn_control%lu_bound_states(2, ispin)
     305              : 
     306          558 :             IF (PRESENT(loc_coeff)) THEN
     307          380 :                mat_ptr = loc_coeff(ispin)
     308              :             ELSE
     309          178 :                mat_ptr => mo_coeff
     310              :             END IF
     311          558 :             IF ((localized_wfn_control%set_of_states == state_loc_list) .OR. &
     312              :                 (localized_wfn_control%set_of_states == state_loc_mixed)) THEN
     313          408 :                ALLOCATE (vecbuffer(1, nao))
     314          136 :                IF (localized_wfn_control%do_homo) THEN
     315          122 :                   my_occ = occupations(localized_wfn_control%loc_states(1, ispin))
     316              :                END IF
     317          136 :                nmosub = SIZE(localized_wfn_control%loc_states, 1)
     318          136 :                CPASSERT(nmosub > 0)
     319          136 :                imoloc = 0
     320          908 :                DO i = lb, ub
     321              :                   ! Get the index in the subset
     322          772 :                   imoloc = imoloc + 1
     323              :                   ! Get the index in the full set
     324          772 :                   imo = localized_wfn_control%loc_states(i, ispin)
     325          772 :                   IF (localized_wfn_control%do_homo) THEN
     326          638 :                      occ_imo = occupations(imo)
     327          638 :                      IF (ABS(occ_imo - my_occ) > localized_wfn_control%eps_occ) THEN
     328            0 :                         IF (localized_wfn_control%localization_method /= do_loc_none) &
     329              :                            CALL cp_abort(__LOCATION__, &
     330              :                                          "States with different occupations "// &
     331            0 :                                          "cannot be rotated together")
     332              :                      END IF
     333              :                   END IF
     334              :                   ! Take the imo vector from the full set and copy in the imoloc vector of the subset
     335              :                   CALL cp_fm_get_submatrix(mat_ptr, vecbuffer, 1, imo, &
     336          772 :                                            nao, 1, transpose=.TRUE.)
     337              :                   CALL cp_fm_set_submatrix(moloc_coeff(ispin), vecbuffer, 1, imoloc, &
     338          908 :                                            nao, 1, transpose=.TRUE.)
     339              :                END DO
     340          136 :                DEALLOCATE (vecbuffer)
     341              :             ELSE
     342          422 :                my_occ = occupations(lb)
     343          422 :                occ_imo = occupations(ub)
     344          422 :                IF (ABS(occ_imo - my_occ) > localized_wfn_control%eps_occ) THEN
     345            0 :                   IF (localized_wfn_control%localization_method /= do_loc_none) &
     346              :                      CALL cp_abort(__LOCATION__, &
     347              :                                    "States with different occupations "// &
     348            0 :                                    "cannot be rotated together")
     349              :                END IF
     350          422 :                nmosub = localized_wfn_control%nloc_states(ispin)
     351          422 :                CALL cp_fm_to_fm(mat_ptr, moloc_coeff(ispin), nmosub, lb, 1)
     352              :             END IF
     353              : 
     354              :             ! we have the mo's to be localized now, see if we can rotate them according to the history
     355              :             ! only do that if we have a history of course. The history is filled
     356         1562 :             IF (PRESENT(mo_loc_history)) THEN
     357          100 :                IF (localized_wfn_control%use_history .AND. ASSOCIATED(mo_loc_history)) THEN
     358              :                   CALL rotate_state_to_ref(moloc_coeff(ispin), &
     359            8 :                                            mo_loc_history(ispin), matrix_s(1)%matrix)
     360              :                END IF
     361              :             END IF
     362              : 
     363              :          END DO
     364              : 
     365          446 :          IF (PRESENT(loc_coeff)) DEALLOCATE (mat_ptr)
     366              : 
     367              :          CALL set_qs_loc_env(qs_loc_env=qs_loc_env, cell=cell, local_molecules=local_molecules, &
     368              :                              moloc_coeff=moloc_coeff, particle_set=particle_set, para_env=para_env, &
     369          446 :                              localized_wfn_control=localized_wfn_control)
     370              : 
     371              :          ! Prepare the operators
     372          446 :          NULLIFY (tmp_fm_struct, mo_coeff)
     373         1338 :          nmosub = MAXVAL(localized_wfn_control%nloc_states)
     374          446 :          CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
     375              :          CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmosub, &
     376          446 :                                   ncol_global=nmosub, para_env=para_env, context=mo_coeff%matrix_struct%context)
     377              : 
     378          446 :          IF (localized_wfn_control%operator_type == op_loc_berry) THEN
     379          446 :             IF (qs_loc_env%cell%orthorhombic) THEN
     380          434 :                dim_op = 3
     381              :             ELSE
     382           12 :                dim_op = 6
     383              :             END IF
     384          446 :             CALL set_qs_loc_env(qs_loc_env=qs_loc_env, dim_op=dim_op)
     385         5460 :             ALLOCATE (qs_loc_env%op_sm_set(2, dim_op))
     386         1820 :             DO i = 1, dim_op
     387         4568 :                DO j = 1, SIZE(qs_loc_env%op_sm_set, 1)
     388         2748 :                   NULLIFY (qs_loc_env%op_sm_set(j, i)%matrix)
     389         2748 :                   ALLOCATE (qs_loc_env%op_sm_set(j, i)%matrix)
     390              :                   CALL dbcsr_copy(qs_loc_env%op_sm_set(j, i)%matrix, matrix_s(1)%matrix, &
     391         2748 :                                   name="qs_loc_env%op_sm_"//TRIM(ADJUSTL(cp_to_string(j)))//"-"//TRIM(ADJUSTL(cp_to_string(i))))
     392         4122 :                   CALL dbcsr_set(qs_loc_env%op_sm_set(j, i)%matrix, 0.0_dp)
     393              :                END DO
     394              :             END DO
     395              : 
     396            0 :          ELSEIF (localized_wfn_control%operator_type == op_loc_pipek) THEN
     397            0 :             natoms = SIZE(qs_loc_env%particle_set, 1)
     398            0 :             ALLOCATE (qs_loc_env%op_fm_set(natoms, 1))
     399            0 :             CALL set_qs_loc_env(qs_loc_env=qs_loc_env, dim_op=natoms)
     400            0 :             DO ispin = 1, SIZE(qs_loc_env%op_fm_set, 2)
     401            0 :                CALL get_mo_set(mos(ispin), nmo=nmo)
     402            0 :                DO iatom = 1, natoms
     403            0 :                   CALL cp_fm_create(qs_loc_env%op_fm_set(iatom, ispin), tmp_fm_struct)
     404              : 
     405            0 :                   CALL cp_fm_get_info(qs_loc_env%op_fm_set(iatom, ispin), nrow_global=nmosub)
     406            0 :                   CPASSERT(nmo >= nmosub)
     407            0 :                   CALL cp_fm_set_all(qs_loc_env%op_fm_set(iatom, ispin), 0.0_dp)
     408              :                END DO ! iatom
     409              :             END DO ! ispin
     410              :          ELSE
     411            0 :             CPABORT("Type of operator not implemented")
     412              :          END IF
     413          446 :          CALL cp_fm_struct_release(tmp_fm_struct)
     414              : 
     415          446 :          IF (localized_wfn_control%operator_type == op_loc_berry) THEN
     416              : 
     417          446 :             CALL initialize_weights(qs_loc_env%cell, qs_loc_env%weights)
     418              : 
     419          446 :             CALL get_berry_operator(qs_loc_env, qs_env)
     420              : 
     421              :          ELSEIF (localized_wfn_control%operator_type == op_loc_pipek) THEN
     422              : 
     423              :             !!    here we don't have to do anything
     424              :             !!    CALL get_pipek_mezey_operator ( qs_loc_env, qs_env )
     425              : 
     426              :          END IF
     427              : 
     428          446 :          qs_loc_env%molecular_states = .FALSE.
     429          446 :          qs_loc_env%wannier_states = .FALSE.
     430              :       END IF
     431          446 :       CALL timestop(handle)
     432              : 
     433          446 :    END SUBROUTINE qs_loc_env_init
     434              : 
     435              : ! **************************************************************************************************
     436              : !> \brief A wrapper to compute the Berry operator for periodic systems
     437              : !> \param qs_loc_env new environment for the localization calculations
     438              : !> \param qs_env the qs_env in which the qs_env lives
     439              : !> \par History
     440              : !>      04.2005 created [MI]
     441              : !>      04.2018 modified [RZK, ZL]
     442              : !> \author MI
     443              : ! **************************************************************************************************
     444          446 :    SUBROUTINE get_berry_operator(qs_loc_env, qs_env)
     445              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
     446              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     447              : 
     448              :       CHARACTER(len=*), PARAMETER :: routineN = 'get_berry_operator'
     449              : 
     450              :       INTEGER                                            :: dim_op, handle
     451              :       TYPE(cell_type), POINTER                           :: cell
     452          446 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: op_sm_set
     453              : 
     454          446 :       CALL timeset(routineN, handle)
     455              : 
     456          446 :       NULLIFY (cell, op_sm_set)
     457              :       CALL get_qs_loc_env(qs_loc_env=qs_loc_env, op_sm_set=op_sm_set, &
     458          446 :                           cell=cell, dim_op=dim_op)
     459          446 :       CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
     460              : 
     461          446 :       CALL timestop(handle)
     462          446 :    END SUBROUTINE get_berry_operator
     463              : 
     464              : ! **************************************************************************************************
     465              : !> \brief Computes the Berry operator for periodic systems
     466              : !>       used to define the spread of the MOS
     467              : !>       Here the matrix elements of the type <mu|cos(kr)|nu> and  <mu|sin(kr)|nu>
     468              : !>       are computed, where mu and nu are the contracted basis functions.
     469              : !>       Namely the Berry operator is exp(ikr)
     470              : !>       k is defined somewhere
     471              : !>       the pair lists are exploited and sparse matrixes are constructed
     472              : !> \param qs_env the qs_env in which the qs_env lives
     473              : !> \param cell ...
     474              : !> \param op_sm_set ...
     475              : !> \param dim_op ...
     476              : !> \par History
     477              : !>      04.2005 created [MI]
     478              : !>      04.2018 wrapped old code [RZK, ZL]
     479              : !> \author MI
     480              : !> \note
     481              : !>      The intgrals are computed analytically  using the primitives GTO
     482              : !>      The contraction is performed block-wise
     483              : ! **************************************************************************************************
     484          472 :    SUBROUTINE compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
     485              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     486              :       TYPE(cell_type), POINTER                           :: cell
     487              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: op_sm_set
     488              :       INTEGER                                            :: dim_op
     489              : 
     490              :       CHARACTER(len=*), PARAMETER :: routineN = 'compute_berry_operator'
     491              : 
     492              :       INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
     493              :          ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nrow, nseta, nsetb, sgfa, sgfb
     494              :       INTEGER, DIMENSION(3)                              :: perd0
     495          472 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     496          472 :                                                             npgfb, nsgfa, nsgfb
     497          472 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     498              :       LOGICAL                                            :: found, new_atom_b
     499              :       REAL(KIND=dp)                                      :: dab, kvec(3), rab2, vector_k(3, 6)
     500              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rb
     501          472 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     502          472 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cosab, rpgfa, rpgfb, sinab, sphi_a, &
     503          472 :                                                             sphi_b, work, zeta, zetb
     504          472 :       TYPE(block_p_type), DIMENSION(:), POINTER          :: op_cos, op_sin
     505          472 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     506              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     507              :       TYPE(neighbor_list_iterator_p_type), &
     508          472 :          DIMENSION(:), POINTER                           :: nl_iterator
     509              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     510          472 :          POINTER                                         :: sab_orb
     511          472 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     512          472 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     513              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     514              : 
     515          472 :       CALL timeset(routineN, handle)
     516          472 :       NULLIFY (qs_kind, qs_kind_set)
     517          472 :       NULLIFY (particle_set)
     518          472 :       NULLIFY (sab_orb)
     519          472 :       NULLIFY (cosab, sinab, work)
     520          472 :       NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
     521          472 :       NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
     522              : 
     523              :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
     524          472 :                       particle_set=particle_set, sab_orb=sab_orb)
     525              : 
     526          472 :       nkind = SIZE(qs_kind_set)
     527              : 
     528              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     529          472 :                            maxco=ldwork, maxlgto=maxl)
     530          472 :       ldab = ldwork
     531         1888 :       ALLOCATE (cosab(ldab, ldab))
     532       128084 :       cosab = 0.0_dp
     533         1416 :       ALLOCATE (sinab(ldab, ldab))
     534       128084 :       sinab = 0.0_dp
     535         1416 :       ALLOCATE (work(ldwork, ldwork))
     536       128084 :       work = 0.0_dp
     537              : 
     538         2868 :       ALLOCATE (op_cos(dim_op))
     539         2396 :       ALLOCATE (op_sin(dim_op))
     540         1924 :       DO i = 1, dim_op
     541         1452 :          NULLIFY (op_cos(i)%block)
     542         1924 :          NULLIFY (op_sin(i)%block)
     543              :       END DO
     544              : 
     545          472 :       kvec = 0.0_dp
     546          472 :       vector_k = 0.0_dp
     547         1888 :       vector_k(:, 1) = twopi*cell%h_inv(1, :)
     548         1888 :       vector_k(:, 2) = twopi*cell%h_inv(2, :)
     549         1888 :       vector_k(:, 3) = twopi*cell%h_inv(3, :)
     550         1888 :       vector_k(:, 4) = twopi*(cell%h_inv(1, :) + cell%h_inv(2, :))
     551         1888 :       vector_k(:, 5) = twopi*(cell%h_inv(1, :) + cell%h_inv(3, :))
     552         1888 :       vector_k(:, 6) = twopi*(cell%h_inv(2, :) + cell%h_inv(3, :))
     553              : 
     554              :       ! This operator can be used only for periodic systems
     555              :       ! If an isolated system is used the periodicity is overimposed
     556         1888 :       perd0(1:3) = cell%perd(1:3)
     557         1888 :       cell%perd(1:3) = 1
     558              : 
     559         2252 :       ALLOCATE (basis_set_list(nkind))
     560         1308 :       DO ikind = 1, nkind
     561          836 :          qs_kind => qs_kind_set(ikind)
     562          836 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
     563         1308 :          IF (ASSOCIATED(basis_set_a)) THEN
     564          836 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     565              :          ELSE
     566            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     567              :          END IF
     568              :       END DO
     569          472 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     570        73599 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     571              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
     572        73127 :                                 iatom=iatom, jatom=jatom, r=rab)
     573        73127 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     574        73127 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     575        73127 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     576        73127 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     577        73127 :          ra = pbc(particle_set(iatom)%r, cell)
     578              :          ! basis ikind
     579        73127 :          first_sgfa => basis_set_a%first_sgf
     580        73127 :          la_max => basis_set_a%lmax
     581        73127 :          la_min => basis_set_a%lmin
     582        73127 :          npgfa => basis_set_a%npgf
     583        73127 :          nseta = basis_set_a%nset
     584        73127 :          nsgfa => basis_set_a%nsgf_set
     585        73127 :          rpgfa => basis_set_a%pgf_radius
     586        73127 :          set_radius_a => basis_set_a%set_radius
     587        73127 :          sphi_a => basis_set_a%sphi
     588        73127 :          zeta => basis_set_a%zet
     589              :          ! basis jkind
     590        73127 :          first_sgfb => basis_set_b%first_sgf
     591        73127 :          lb_max => basis_set_b%lmax
     592        73127 :          lb_min => basis_set_b%lmin
     593        73127 :          npgfb => basis_set_b%npgf
     594        73127 :          nsetb = basis_set_b%nset
     595        73127 :          nsgfb => basis_set_b%nsgf_set
     596        73127 :          rpgfb => basis_set_b%pgf_radius
     597        73127 :          set_radius_b => basis_set_b%set_radius
     598        73127 :          sphi_b => basis_set_b%sphi
     599        73127 :          zetb => basis_set_b%zet
     600              : 
     601        73127 :          ldsa = SIZE(sphi_a, 1)
     602        73127 :          ldsb = SIZE(sphi_b, 1)
     603        73127 :          IF (inode == 1) last_jatom = 0
     604              : 
     605       292508 :          rb = rab + ra
     606              : 
     607        73127 :          IF (jatom /= last_jatom) THEN
     608              :             new_atom_b = .TRUE.
     609              :             last_jatom = jatom
     610              :          ELSE
     611              :             new_atom_b = .FALSE.
     612              :          END IF
     613              : 
     614              :          IF (new_atom_b) THEN
     615        17801 :             IF (iatom <= jatom) THEN
     616         9298 :                irow = iatom
     617         9298 :                icol = jatom
     618              :             ELSE
     619         8503 :                irow = jatom
     620         8503 :                icol = iatom
     621              :             END IF
     622              : 
     623        71672 :             DO i = 1, dim_op
     624        53871 :                NULLIFY (op_cos(i)%block)
     625              :                CALL dbcsr_get_block_p(matrix=op_sm_set(1, i)%matrix, &
     626        53871 :                                       row=irow, col=icol, block=op_cos(i)%block, found=found)
     627        53871 :                NULLIFY (op_sin(i)%block)
     628              :                CALL dbcsr_get_block_p(matrix=op_sm_set(2, i)%matrix, &
     629        71672 :                                       row=irow, col=icol, block=op_sin(i)%block, found=found)
     630              :             END DO
     631              :          END IF ! new_atom_b
     632              : 
     633        73127 :          rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
     634        73127 :          dab = SQRT(rab2)
     635              : 
     636        73127 :          nrow = 0
     637       229459 :          DO iset = 1, nseta
     638              : 
     639       155860 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     640       155860 :             sgfa = first_sgfa(1, iset)
     641              : 
     642       508274 :             DO jset = 1, nsetb
     643              : 
     644       352414 :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     645       352414 :                sgfb = first_sgfb(1, jset)
     646              : 
     647       508274 :                IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
     648              : 
     649              : !           *** Calculate the primitive overlap integrals ***
     650       717282 :                   DO i = 1, dim_op
     651      2223564 :                      kvec(1:3) = vector_k(1:3, i)
     652    188199543 :                      cosab = 0.0_dp
     653    188199543 :                      sinab = 0.0_dp
     654              :                      CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
     655              :                                  la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
     656              :                                  rpgfb(:, jset), lb_min(jset), &
     657       555891 :                                  ra, rb, kvec, cosab, sinab)
     658              :                      CALL contract_cossin(op_cos(i)%block, op_sin(i)%block, &
     659              :                                           iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
     660              :                                           jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
     661       717282 :                                           cosab, sinab, ldab, work, ldwork)
     662              :                   END DO
     663              : 
     664              :                END IF !  >= dab
     665              : 
     666              :             END DO ! jset
     667              : 
     668       228987 :             nrow = nrow + ncoa
     669              : 
     670              :          END DO ! iset
     671              : 
     672              :       END DO
     673          472 :       CALL neighbor_list_iterator_release(nl_iterator)
     674              : 
     675              :       ! Set back the correct periodicity
     676         1888 :       cell%perd(1:3) = perd0(1:3)
     677              : 
     678         1924 :       DO i = 1, dim_op
     679         1452 :          NULLIFY (op_cos(i)%block)
     680         1924 :          NULLIFY (op_sin(i)%block)
     681              :       END DO
     682          472 :       DEALLOCATE (op_cos, op_sin)
     683              : 
     684          472 :       DEALLOCATE (cosab, sinab, work, basis_set_list)
     685              : 
     686          472 :       CALL timestop(handle)
     687          944 :    END SUBROUTINE compute_berry_operator
     688              : 
     689              : ! **************************************************************************************************
     690              : !> \brief ...
     691              : !> \param qs_loc_env ...
     692              : !> \param section ...
     693              : !> \param mo_array ...
     694              : !> \param coeff_localized ...
     695              : !> \param do_homo ...
     696              : !> \param evals ...
     697              : !> \param do_mixed ...
     698              : ! **************************************************************************************************
     699          294 :    SUBROUTINE loc_write_restart(qs_loc_env, section, mo_array, coeff_localized, &
     700              :                                 do_homo, evals, do_mixed)
     701              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
     702              :       TYPE(section_vals_type), POINTER                   :: section
     703              :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     704              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: coeff_localized
     705              :       LOGICAL, INTENT(IN)                                :: do_homo
     706              :       TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
     707              :          POINTER                                         :: evals
     708              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_mixed
     709              : 
     710              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'loc_write_restart'
     711              : 
     712              :       CHARACTER(LEN=default_path_length)                 :: filename
     713              :       CHARACTER(LEN=default_string_length)               :: my_middle
     714              :       INTEGER                                            :: handle, ispin, max_block, nao, nloc, &
     715              :                                                             nmo, output_unit, rst_unit
     716              :       LOGICAL                                            :: my_do_mixed
     717              :       TYPE(cp_logger_type), POINTER                      :: logger
     718              :       TYPE(section_vals_type), POINTER                   :: print_key
     719              : 
     720          294 :       CALL timeset(routineN, handle)
     721          294 :       NULLIFY (logger)
     722          294 :       logger => cp_get_default_logger()
     723          294 :       output_unit = cp_logger_get_default_io_unit(logger)
     724              : 
     725          294 :       IF (qs_loc_env%do_localize) THEN
     726              : 
     727          278 :          print_key => section_vals_get_subs_vals(section, "LOC_RESTART")
     728          278 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     729              :                                               section, "LOC_RESTART"), &
     730              :                    cp_p_file)) THEN
     731              : 
     732              :             ! Open file
     733           30 :             rst_unit = -1
     734              : 
     735           30 :             my_do_mixed = .FALSE.
     736           30 :             IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
     737           30 :             IF (do_homo) THEN
     738           30 :                my_middle = "LOC_HOMO"
     739            0 :             ELSEIF (my_do_mixed) THEN
     740            0 :                my_middle = "LOC_MIXED"
     741              :             ELSE
     742            0 :                my_middle = "LOC_LUMO"
     743              :             END IF
     744              : 
     745              :             rst_unit = cp_print_key_unit_nr(logger, section, "LOC_RESTART", &
     746              :                                             extension=".wfn", file_status="REPLACE", file_action="WRITE", &
     747           30 :                                             file_form="UNFORMATTED", middle_name=TRIM(my_middle))
     748              : 
     749           30 :             IF (rst_unit > 0) filename = cp_print_key_generate_filename(logger, print_key, &
     750              :                                                                         middle_name=TRIM(my_middle), extension=".wfn", &
     751           15 :                                                                         my_local=.FALSE.)
     752              : 
     753           30 :             IF (output_unit > 0) THEN
     754              :                WRITE (UNIT=output_unit, FMT="(/,T2,A, A/)") &
     755           15 :                   "LOCALIZATION| Write restart file for the localized MOS : ", &
     756           30 :                   TRIM(filename)
     757              :             END IF
     758              : 
     759           30 :             IF (rst_unit > 0) THEN
     760           15 :                WRITE (rst_unit) qs_loc_env%localized_wfn_control%set_of_states
     761           15 :                WRITE (rst_unit) qs_loc_env%localized_wfn_control%lu_bound_states
     762           15 :                WRITE (rst_unit) qs_loc_env%localized_wfn_control%nloc_states
     763              :             END IF
     764              : 
     765           70 :             DO ispin = 1, SIZE(coeff_localized)
     766           30 :                ASSOCIATE (mo_coeff => coeff_localized(ispin))
     767           40 :                   CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_global=nmo, ncol_block=max_block)
     768           40 :                   nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
     769           40 :                   IF (rst_unit > 0) THEN
     770           20 :                      WRITE (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
     771           20 :                      IF (do_homo .OR. my_do_mixed) THEN
     772           20 :                         WRITE (rst_unit) nmo, &
     773           20 :                            mo_array(ispin)%homo, &
     774           20 :                            mo_array(ispin)%lfomo, &
     775           40 :                            mo_array(ispin)%nelectron
     776          456 :                         WRITE (rst_unit) mo_array(ispin)%eigenvalues(1:nmo), &
     777          476 :                            mo_array(ispin)%occupation_numbers(1:nmo)
     778              :                      ELSE
     779            0 :                         WRITE (rst_unit) nmo
     780            0 :                         WRITE (rst_unit) evals(ispin)%array(1:nmo)
     781              :                      END IF
     782              :                   END IF
     783              : 
     784           80 :                   CALL cp_fm_write_unformatted(mo_coeff, rst_unit)
     785              :                END ASSOCIATE
     786              : 
     787              :             END DO
     788              : 
     789              :             CALL cp_print_key_finished_output(rst_unit, logger, section, &
     790           30 :                                               "LOC_RESTART")
     791              :          END IF
     792              : 
     793              :       END IF
     794              : 
     795          294 :       CALL timestop(handle)
     796              : 
     797          294 :    END SUBROUTINE loc_write_restart
     798              : 
     799              : ! **************************************************************************************************
     800              : !> \brief ...
     801              : !> \param qs_loc_env ...
     802              : !> \param mos ...
     803              : !> \param mos_localized ...
     804              : !> \param section ...
     805              : !> \param section2 ...
     806              : !> \param para_env ...
     807              : !> \param do_homo ...
     808              : !> \param restart_found ...
     809              : !> \param evals ...
     810              : !> \param do_mixed ...
     811              : ! **************************************************************************************************
     812            6 :    SUBROUTINE loc_read_restart(qs_loc_env, mos, mos_localized, section, section2, para_env, &
     813              :                                do_homo, restart_found, evals, do_mixed)
     814              : 
     815              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
     816              :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     817              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: mos_localized
     818              :       TYPE(section_vals_type), POINTER                   :: section, section2
     819              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     820              :       LOGICAL, INTENT(IN)                                :: do_homo
     821              :       LOGICAL, INTENT(INOUT)                             :: restart_found
     822              :       TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
     823              :          POINTER                                         :: evals
     824              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_mixed
     825              : 
     826              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'loc_read_restart'
     827              : 
     828              :       CHARACTER(LEN=25)                                  :: fname_key
     829              :       CHARACTER(LEN=default_path_length)                 :: filename
     830              :       CHARACTER(LEN=default_string_length)               :: my_middle
     831              :       INTEGER :: handle, homo_read, i, ispin, lfomo_read, max_nloc, n_rep_val, nao, &
     832              :          nelectron_read, nloc, nmo, nmo_read, nspin, output_unit, rst_unit
     833              :       LOGICAL                                            :: file_exists, my_do_mixed
     834            6 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eig_read, occ_read
     835            6 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
     836              :       TYPE(cp_logger_type), POINTER                      :: logger
     837              :       TYPE(section_vals_type), POINTER                   :: print_key
     838              : 
     839            6 :       CALL timeset(routineN, handle)
     840              : 
     841            6 :       logger => cp_get_default_logger()
     842              : 
     843            6 :       nspin = SIZE(mos_localized)
     844            6 :       nao = mos(1)%nao
     845            6 :       rst_unit = -1
     846              : 
     847              :       output_unit = cp_print_key_unit_nr(logger, section2, &
     848            6 :                                          "PROGRAM_RUN_INFO", extension=".Log")
     849              : 
     850            6 :       my_do_mixed = .FALSE.
     851            6 :       IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
     852            6 :       IF (do_homo) THEN
     853            6 :          fname_key = "LOCHOMO_RESTART_FILE_NAME"
     854            0 :       ELSEIF (my_do_mixed) THEN
     855            0 :          fname_key = "LOCMIXD_RESTART_FILE_NAME"
     856              :       ELSE
     857            0 :          fname_key = "LOCLUMO_RESTART_FILE_NAME"
     858            0 :          IF (.NOT. PRESENT(evals)) &
     859            0 :             CPABORT("Missing argument to localize unoccupied states.")
     860              :       END IF
     861              : 
     862            6 :       file_exists = .FALSE.
     863            6 :       CALL section_vals_val_get(section, fname_key, n_rep_val=n_rep_val)
     864            6 :       IF (n_rep_val > 0) THEN
     865            0 :          CALL section_vals_val_get(section, fname_key, c_val=filename)
     866              :       ELSE
     867              : 
     868            6 :          print_key => section_vals_get_subs_vals(section2, "LOC_RESTART")
     869            6 :          IF (do_homo) THEN
     870            6 :             my_middle = "LOC_HOMO"
     871            0 :          ELSEIF (my_do_mixed) THEN
     872            0 :             my_middle = "LOC_MIXED"
     873              :          ELSE
     874            0 :             my_middle = "LOC_LUMO"
     875              :          END IF
     876              :          filename = cp_print_key_generate_filename(logger, print_key, &
     877              :                                                    middle_name=TRIM(my_middle), extension=".wfn", &
     878            6 :                                                    my_local=.FALSE.)
     879              :       END IF
     880              : 
     881            6 :       IF (para_env%is_source()) INQUIRE (FILE=filename, exist=file_exists)
     882              : 
     883            6 :       IF (file_exists) THEN
     884            2 :          IF (para_env%is_source()) THEN
     885              :             CALL open_file(file_name=filename, &
     886              :                            file_action="READ", &
     887              :                            file_form="UNFORMATTED", &
     888              :                            file_status="OLD", &
     889            2 :                            unit_number=rst_unit)
     890              : 
     891            2 :             READ (rst_unit) qs_loc_env%localized_wfn_control%set_of_states
     892            2 :             READ (rst_unit) qs_loc_env%localized_wfn_control%lu_bound_states
     893            2 :             READ (rst_unit) qs_loc_env%localized_wfn_control%nloc_states
     894              :          END IF
     895              :       ELSE
     896            4 :          IF (output_unit > 0) WRITE (output_unit, "(/,T10,A)") &
     897            1 :             "Restart file not available filename=<"//TRIM(filename)//'>'
     898              :       END IF
     899            6 :       CALL para_env%bcast(file_exists)
     900              : 
     901            6 :       IF (file_exists) THEN
     902            4 :          restart_found = .TRUE.
     903              : 
     904            4 :          CALL para_env%bcast(qs_loc_env%localized_wfn_control%set_of_states)
     905            4 :          CALL para_env%bcast(qs_loc_env%localized_wfn_control%lu_bound_states)
     906            4 :          CALL para_env%bcast(qs_loc_env%localized_wfn_control%nloc_states)
     907              : 
     908           12 :          max_nloc = MAXVAL(qs_loc_env%localized_wfn_control%nloc_states(:))
     909              : 
     910           12 :          ALLOCATE (vecbuffer(1, nao))
     911            4 :          IF (ASSOCIATED(qs_loc_env%localized_wfn_control%loc_states)) THEN
     912            2 :             DEALLOCATE (qs_loc_env%localized_wfn_control%loc_states)
     913              :          END IF
     914           12 :          ALLOCATE (qs_loc_env%localized_wfn_control%loc_states(max_nloc, 2))
     915           56 :          qs_loc_env%localized_wfn_control%loc_states = 0
     916              : 
     917           10 :          DO ispin = 1, nspin
     918            6 :             IF (do_homo .OR. do_mixed) THEN
     919            6 :                nmo = mos(ispin)%nmo
     920              :             ELSE
     921            0 :                nmo = SIZE(evals(ispin)%array, 1)
     922              :             END IF
     923            6 :             IF (para_env%is_source() .AND. (nmo > 0)) THEN
     924            3 :                nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
     925            3 :                READ (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
     926            3 :                IF (do_homo .OR. do_mixed) THEN
     927            3 :                   READ (rst_unit) nmo_read, homo_read, lfomo_read, nelectron_read
     928           12 :                   ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
     929           79 :                   eig_read = 0.0_dp
     930           79 :                   occ_read = 0.0_dp
     931            3 :                   READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read)
     932              :                ELSE
     933            0 :                   READ (rst_unit) nmo_read
     934            0 :                   ALLOCATE (eig_read(nmo_read))
     935            0 :                   eig_read = 0.0_dp
     936            0 :                   READ (rst_unit) eig_read(1:nmo_read)
     937              :                END IF
     938            3 :                IF (nmo_read < nmo) &
     939              :                   CALL cp_warn(__LOCATION__, &
     940              :                                "The number of MOs on the restart unit is smaller than the number of "// &
     941            0 :                                "the allocated MOs. ")
     942            3 :                IF (nmo_read > nmo) &
     943              :                   CALL cp_warn(__LOCATION__, &
     944              :                                "The number of MOs on the restart unit is greater than the number of "// &
     945            0 :                                "the allocated MOs. The read MO set will be truncated!")
     946              : 
     947            3 :                nmo = MIN(nmo, nmo_read)
     948            3 :                IF (do_homo .OR. do_mixed) THEN
     949           79 :                   mos(ispin)%eigenvalues(1:nmo) = eig_read(1:nmo)
     950           79 :                   mos(ispin)%occupation_numbers(1:nmo) = occ_read(1:nmo)
     951            3 :                   DEALLOCATE (eig_read, occ_read)
     952              :                ELSE
     953            0 :                   evals(ispin)%array(1:nmo) = eig_read(1:nmo)
     954            0 :                   DEALLOCATE (eig_read)
     955              :                END IF
     956              : 
     957              :             END IF
     958            6 :             IF (do_homo .OR. do_mixed) THEN
     959          310 :                CALL para_env%bcast(mos(ispin)%eigenvalues)
     960          310 :                CALL para_env%bcast(mos(ispin)%occupation_numbers)
     961              :             ELSE
     962            0 :                CALL para_env%bcast(evals(ispin)%array)
     963              :             END IF
     964              : 
     965          162 :             DO i = 1, nmo
     966          152 :                IF (para_env%is_source()) THEN
     967        15236 :                   READ (rst_unit) vecbuffer
     968              :                ELSE
     969         7656 :                   vecbuffer(1, :) = 0.0_dp
     970              :                END IF
     971        60792 :                CALL para_env%bcast(vecbuffer)
     972              :                CALL cp_fm_set_submatrix(mos_localized(ispin), &
     973          158 :                                         vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
     974              :             END DO
     975              :          END DO
     976              : 
     977          108 :          CALL para_env%bcast(qs_loc_env%localized_wfn_control%loc_states)
     978              : 
     979            4 :          DEALLOCATE (vecbuffer)
     980              : 
     981              :       END IF
     982              : 
     983              :       ! Close restart file
     984            6 :       IF (para_env%is_source()) THEN
     985            3 :          IF (file_exists) CALL close_file(unit_number=rst_unit)
     986              :       END IF
     987              : 
     988            6 :       CALL timestop(handle)
     989              : 
     990           12 :    END SUBROUTINE loc_read_restart
     991              : 
     992              : ! **************************************************************************************************
     993              : !> \brief initializes everything needed for localization of the HOMOs
     994              : !> \param qs_loc_env ...
     995              : !> \param loc_section ...
     996              : !> \param do_homo ...
     997              : !> \param do_mixed ...
     998              : !> \param do_xas ...
     999              : !> \param nloc_xas ...
    1000              : !> \param spin_xas ...
    1001              : !> \par History
    1002              : !>      2009 created
    1003              : ! **************************************************************************************************
    1004          394 :    SUBROUTINE qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, &
    1005              :                                   do_xas, nloc_xas, spin_xas)
    1006              : 
    1007              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
    1008              :       TYPE(section_vals_type), POINTER                   :: loc_section
    1009              :       LOGICAL, INTENT(IN)                                :: do_homo
    1010              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_mixed, do_xas
    1011              :       INTEGER, INTENT(IN), OPTIONAL                      :: nloc_xas, spin_xas
    1012              : 
    1013              :       LOGICAL                                            :: my_do_mixed
    1014              :       TYPE(localized_wfn_control_type), POINTER          :: localized_wfn_control
    1015              : 
    1016          394 :       NULLIFY (localized_wfn_control)
    1017              : 
    1018          394 :       IF (PRESENT(do_mixed)) THEN
    1019            2 :          my_do_mixed = do_mixed
    1020              :       ELSE
    1021          392 :          my_do_mixed = .FALSE.
    1022              :       END IF
    1023          394 :       CALL localized_wfn_control_create(localized_wfn_control)
    1024          394 :       CALL set_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
    1025          394 :       CALL localized_wfn_control_release(localized_wfn_control)
    1026          394 :       CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
    1027          394 :       localized_wfn_control%do_homo = do_homo
    1028          394 :       localized_wfn_control%do_mixed = my_do_mixed
    1029              :       CALL read_loc_section(localized_wfn_control, loc_section, qs_loc_env%do_localize, &
    1030          394 :                             my_do_mixed, do_xas, nloc_xas, spin_xas)
    1031              : 
    1032          394 :    END SUBROUTINE qs_loc_control_init
    1033              : 
    1034              : ! **************************************************************************************************
    1035              : !> \brief initializes everything needed for localization of the molecular orbitals
    1036              : !> \param qs_env ...
    1037              : !> \param qs_loc_env ...
    1038              : !> \param localize_section ...
    1039              : !> \param mos_localized ...
    1040              : !> \param do_homo ...
    1041              : !> \param do_mo_cubes ...
    1042              : !> \param mo_loc_history ...
    1043              : !> \param evals ...
    1044              : !> \param tot_zeff_corr ...
    1045              : !> \param do_mixed ...
    1046              : ! **************************************************************************************************
    1047          294 :    SUBROUTINE qs_loc_init(qs_env, qs_loc_env, localize_section, mos_localized, &
    1048              :                           do_homo, do_mo_cubes, mo_loc_history, evals, &
    1049              :                           tot_zeff_corr, do_mixed)
    1050              : 
    1051              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1052              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
    1053              :       TYPE(section_vals_type), POINTER                   :: localize_section
    1054              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: mos_localized
    1055              :       LOGICAL, OPTIONAL                                  :: do_homo, do_mo_cubes
    1056              :       TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER  :: mo_loc_history
    1057              :       TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
    1058              :          POINTER                                         :: evals
    1059              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: tot_zeff_corr
    1060              :       LOGICAL, OPTIONAL                                  :: do_mixed
    1061              : 
    1062              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_loc_init'
    1063              : 
    1064              :       INTEGER :: handle, homo, i, ilast_intocc, ilow, ispin, iup, n_mo(2), n_mos(2), nao, &
    1065              :          nelectron, nextra, nmoloc(2), nocc, npocc, nspin, output_unit
    1066              :       LOGICAL                                            :: my_do_homo, my_do_mixed, my_do_mo_cubes, &
    1067              :                                                             restart_found
    1068              :       REAL(KIND=dp)                                      :: maxocc, my_tot_zeff_corr
    1069          294 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues, occupation
    1070              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1071              :       TYPE(cp_logger_type), POINTER                      :: logger
    1072          294 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_rmpv, mo_derivs
    1073              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1074              :       TYPE(localized_wfn_control_type), POINTER          :: localized_wfn_control
    1075          294 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1076              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1077              :       TYPE(scf_control_type), POINTER                    :: scf_control
    1078              :       TYPE(section_vals_type), POINTER                   :: loc_print_section
    1079              : 
    1080          294 :       CALL timeset(routineN, handle)
    1081              : 
    1082          294 :       NULLIFY (mos, mo_coeff, mo_eigenvalues, occupation, ks_rmpv, mo_derivs, scf_control, para_env)
    1083              :       CALL get_qs_env(qs_env, &
    1084              :                       mos=mos, &
    1085              :                       matrix_ks=ks_rmpv, &
    1086              :                       mo_derivs=mo_derivs, &
    1087              :                       scf_control=scf_control, &
    1088              :                       dft_control=dft_control, &
    1089          294 :                       para_env=para_env)
    1090              : 
    1091          294 :       loc_print_section => section_vals_get_subs_vals(localize_section, "PRINT")
    1092              : 
    1093          294 :       logger => cp_get_default_logger()
    1094          294 :       output_unit = cp_logger_get_default_io_unit(logger)
    1095              : 
    1096          294 :       nspin = SIZE(mos)
    1097          294 :       IF (PRESENT(do_homo)) THEN
    1098          294 :          my_do_homo = do_homo
    1099              :       ELSE
    1100            0 :          my_do_homo = .TRUE.
    1101              :       END IF
    1102          294 :       IF (PRESENT(do_mo_cubes)) THEN
    1103          104 :          my_do_mo_cubes = do_mo_cubes
    1104              :       ELSE
    1105              :          my_do_mo_cubes = .FALSE.
    1106              :       END IF
    1107          294 :       IF (PRESENT(do_mixed)) THEN
    1108            2 :          my_do_mixed = do_mixed
    1109              :       ELSE
    1110          292 :          my_do_mixed = .FALSE.
    1111              :       END IF
    1112          294 :       IF (PRESENT(tot_zeff_corr)) THEN
    1113            2 :          my_tot_zeff_corr = tot_zeff_corr
    1114              :       ELSE
    1115          292 :          my_tot_zeff_corr = 0.0_dp
    1116              :       END IF
    1117          294 :       restart_found = .FALSE.
    1118              : 
    1119          294 :       IF (qs_loc_env%do_localize) THEN
    1120              :          ! Some setup for MOs to be localized
    1121          278 :          CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
    1122          278 :          IF (localized_wfn_control%loc_restart) THEN
    1123            6 :             IF (localized_wfn_control%nextra > 0) THEN
    1124              :                ! currently only the occupied guess is read
    1125            0 :                my_do_homo = .FALSE.
    1126              :             END IF
    1127              :             CALL loc_read_restart(qs_loc_env, mos, mos_localized, localize_section, &
    1128              :                                   loc_print_section, para_env, my_do_homo, restart_found, evals=evals, &
    1129            6 :                                   do_mixed=my_do_mixed)
    1130            9 :             IF (output_unit > 0) WRITE (output_unit, "(/,T2,A,A)") "LOCALIZATION| ", &
    1131            6 :                "   The orbitals to be localized are read from localization restart file."
    1132           18 :             nmoloc = localized_wfn_control%nloc_states
    1133           18 :             localized_wfn_control%nguess = nmoloc
    1134            6 :             IF (localized_wfn_control%nextra > 0) THEN
    1135              :                ! reset different variables in localized_wfn_control:
    1136              :                ! lu_bound_states, nloc_states, loc_states
    1137            0 :                localized_wfn_control%loc_restart = restart_found
    1138            0 :                localized_wfn_control%set_of_states = state_loc_mixed
    1139            0 :                DO ispin = 1, nspin
    1140              :                   CALL get_mo_set(mos(ispin), homo=homo, occupation_numbers=occupation, &
    1141            0 :                                   maxocc=maxocc)
    1142            0 :                   nextra = localized_wfn_control%nextra
    1143            0 :                   nocc = homo
    1144            0 :                   DO i = nocc, 1, -1
    1145            0 :                      IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
    1146            0 :                         ilast_intocc = i
    1147            0 :                         EXIT
    1148              :                      END IF
    1149              :                   END DO
    1150            0 :                   nocc = ilast_intocc
    1151            0 :                   npocc = homo - nocc
    1152            0 :                   nmoloc(ispin) = nocc + nextra
    1153            0 :                   localized_wfn_control%lu_bound_states(1, ispin) = 1
    1154            0 :                   localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
    1155            0 :                   localized_wfn_control%nloc_states(ispin) = nmoloc(ispin)
    1156              :                END DO
    1157            0 :                my_do_homo = .FALSE.
    1158              :             END IF
    1159              :          END IF
    1160          278 :          IF (.NOT. restart_found) THEN
    1161          274 :             nmoloc = 0
    1162          648 :             DO ispin = 1, nspin
    1163              :                CALL get_mo_set(mos(ispin), nmo=n_mo(ispin), nelectron=nelectron, homo=homo, nao=nao, &
    1164              :                                mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, occupation_numbers=occupation, &
    1165          374 :                                maxocc=maxocc)
    1166              :                ! Get eigenstates (only needed if not already calculated before)
    1167              :                IF ((.NOT. my_do_mo_cubes) &
    1168              :                    !                  .OR. section_get_ival(dft_section,"PRINT%MO_CUBES%NHOMO")==0)&
    1169          374 :                    .AND. my_do_homo .AND. qs_env%scf_env%method == ot_method_nr .AND. (.NOT. dft_control%restricted)) THEN
    1170           30 :                   CALL make_mo_eig(mos, nspin, ks_rmpv, scf_control, mo_derivs)
    1171              :                END IF
    1172         1022 :                IF (localized_wfn_control%set_of_states == state_loc_all .AND. my_do_homo) THEN
    1173          334 :                   nmoloc(ispin) = NINT(nelectron/occupation(1))
    1174          334 :                   IF (n_mo(ispin) > homo) THEN
    1175           30 :                      DO i = nmoloc(ispin), 1, -1
    1176           30 :                         IF (occupation(1) - occupation(i) < localized_wfn_control%eps_occ) THEN
    1177           14 :                            ilast_intocc = i
    1178           14 :                            EXIT
    1179              :                         END IF
    1180              :                      END DO
    1181              :                   ELSE
    1182          320 :                      ilast_intocc = nmoloc(ispin)
    1183              :                   END IF
    1184          334 :                   nmoloc(ispin) = ilast_intocc
    1185          334 :                   localized_wfn_control%lu_bound_states(1, ispin) = 1
    1186          334 :                   localized_wfn_control%lu_bound_states(2, ispin) = ilast_intocc
    1187          334 :                   IF (nmoloc(ispin) /= n_mo(ispin)) THEN
    1188           14 :                      IF (output_unit > 0) &
    1189              :                         WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,F12.6,A,F12.6,A)") &
    1190            7 :                         "LOCALIZATION| Spin ", ispin, " The first ", &
    1191            7 :                         ilast_intocc, " occupied orbitals are localized,", " with energies from ", &
    1192           14 :                         mo_eigenvalues(1), " to ", mo_eigenvalues(ilast_intocc), " [a.u.]."
    1193              :                   END IF
    1194           40 :                ELSE IF (localized_wfn_control%set_of_states == energy_loc_range .AND. my_do_homo) THEN
    1195           12 :                   ilow = 0
    1196           12 :                   iup = 0
    1197           20 :                   DO i = 1, n_mo(ispin)
    1198           20 :                      IF (mo_eigenvalues(i) >= localized_wfn_control%lu_ene_bound(1)) THEN
    1199              :                         ilow = i
    1200              :                         EXIT
    1201              :                      END IF
    1202              :                   END DO
    1203          324 :                   DO i = n_mo(ispin), 1, -1
    1204          324 :                      IF (mo_eigenvalues(i) <= localized_wfn_control%lu_ene_bound(2)) THEN
    1205              :                         iup = i
    1206              :                         EXIT
    1207              :                      END IF
    1208              :                   END DO
    1209           12 :                   localized_wfn_control%lu_bound_states(1, ispin) = ilow
    1210           12 :                   localized_wfn_control%lu_bound_states(2, ispin) = iup
    1211           12 :                   localized_wfn_control%nloc_states(ispin) = iup - ilow + 1
    1212           12 :                   nmoloc(ispin) = localized_wfn_control%nloc_states(ispin)
    1213           12 :                   IF (occupation(ilow) - occupation(iup) > localized_wfn_control%eps_occ) THEN
    1214              :                      CALL cp_abort(__LOCATION__, &
    1215              :                                    "The selected energy range includes orbitals with different occupation number. "// &
    1216            0 :                                    "The localization procedure cannot be applied.")
    1217              :                   END IF
    1218           18 :                   IF (output_unit > 0) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", ispin, " : ", &
    1219           12 :                      nmoloc(ispin), " orbitals in the selected energy range are localized."
    1220           28 :                ELSE IF (localized_wfn_control%set_of_states == state_loc_all .AND. (.NOT. my_do_homo)) THEN
    1221            0 :                   nmoloc(ispin) = n_mo(ispin) - homo
    1222            0 :                   localized_wfn_control%lu_bound_states(1, ispin) = homo + 1
    1223            0 :                   localized_wfn_control%lu_bound_states(2, ispin) = n_mo(ispin)
    1224            0 :                   IF (output_unit > 0) &
    1225              :                      WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,F12.6,A,F12.6,A)") &
    1226            0 :                      "LOCALIZATION| Spin ", ispin, " The first ", &
    1227            0 :                      nmoloc(ispin), " virtual orbitals are localized,", " with energies from ", &
    1228            0 :                      mo_eigenvalues(homo + 1), " to ", mo_eigenvalues(n_mo(ispin)), " [a.u.]."
    1229           28 :                ELSE IF (localized_wfn_control%set_of_states == state_loc_mixed) THEN
    1230            2 :                   nextra = localized_wfn_control%nextra
    1231            2 :                   nocc = homo
    1232            6 :                   DO i = nocc, 1, -1
    1233            6 :                      IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
    1234            2 :                         ilast_intocc = i
    1235            2 :                         EXIT
    1236              :                      END IF
    1237              :                   END DO
    1238            2 :                   nocc = ilast_intocc
    1239            2 :                   npocc = homo - nocc
    1240            2 :                   nmoloc(ispin) = nocc + nextra
    1241            2 :                   localized_wfn_control%lu_bound_states(1, ispin) = 1
    1242            2 :                   localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
    1243            2 :                   IF (output_unit > 0) &
    1244              :                      WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,I6,/,T15,A,I6,/,T15,A,I6,/,T15,A,F12.6,A)") &
    1245            1 :                      "LOCALIZATION| Spin ", ispin, " The first ", &
    1246            1 :                      nmoloc(ispin), " orbitals are localized.", &
    1247            1 :                      "Number of fully occupied MOs: ", nocc, &
    1248            1 :                      "Number of partially occupied MOs: ", npocc, &
    1249            1 :                      "Number of extra degrees of freedom: ", nextra, &
    1250            2 :                      "Excess charge: ", my_tot_zeff_corr, " electrons"
    1251              :                ELSE
    1252           26 :                   nmoloc(ispin) = MIN(localized_wfn_control%nloc_states(1), n_mo(ispin))
    1253           33 :                   IF (output_unit > 0 .AND. my_do_homo) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", ispin, &
    1254           14 :                      " : ", nmoloc(ispin), " occupied orbitals are localized, as given in the input list."
    1255           19 :                   IF (output_unit > 0 .AND. (.NOT. my_do_homo)) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", &
    1256           12 :                      ispin, " : ", nmoloc(ispin), " unoccupied orbitals are localized, as given in the input list."
    1257           26 :                   IF (n_mo(ispin) > homo .AND. my_do_homo) THEN
    1258            8 :                      ilow = localized_wfn_control%loc_states(1, ispin)
    1259           56 :                      DO i = 2, nmoloc(ispin)
    1260           48 :                         iup = localized_wfn_control%loc_states(i, ispin)
    1261           56 :                         IF (ABS(occupation(ilow) - occupation(iup)) > localized_wfn_control%eps_occ) THEN
    1262              :                            ! write warning
    1263              :                            CALL cp_warn(__LOCATION__, &
    1264              :                                         "User requested the calculation of localized wavefunction from a subset of MOs, "// &
    1265              :                                         "including MOs with different occupations. Check the selected subset, "// &
    1266              :                                         "the electronic density is not invariant with "// &
    1267            0 :                                         "respect to rotations among orbitals with different occupation numbers!")
    1268              :                         END IF
    1269              :                      END DO
    1270              :                   END IF
    1271              :                END IF
    1272              :             END DO ! ispin
    1273          822 :             n_mos(:) = nao - n_mo(:)
    1274          274 :             IF (my_do_homo .OR. my_do_mixed) n_mos = n_mo
    1275          274 :             CALL set_loc_wfn_lists(localized_wfn_control, nmoloc, n_mos, nspin)
    1276              :          END IF
    1277          278 :          CALL set_loc_centers(localized_wfn_control, nmoloc, nspin)
    1278          278 :          IF (my_do_homo .OR. my_do_mixed) THEN
    1279              :             CALL qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, &
    1280          272 :                                  loc_coeff=mos_localized, mo_loc_history=mo_loc_history)
    1281              :          END IF
    1282              :       ELSE
    1283              :          ! Let's inform in case the section is not present in the input
    1284              :          CALL cp_warn(__LOCATION__, &
    1285              :                       "User requested the calculation of the localized wavefunction but the section "// &
    1286           16 :                       "LOCALIZE was not specified. Localization will not be performed!")
    1287              :       END IF
    1288              : 
    1289          294 :       CALL timestop(handle)
    1290              : 
    1291          294 :    END SUBROUTINE qs_loc_init
    1292              : 
    1293              : ! **************************************************************************************************
    1294              : !> \brief read the controlparameter from input, using the new input scheme
    1295              : !> \param localized_wfn_control ...
    1296              : !> \param loc_section ...
    1297              : !> \param localize ...
    1298              : !> \param do_mixed ...
    1299              : !> \param do_xas ...
    1300              : !> \param nloc_xas ...
    1301              : !> \param spin_channel_xas ...
    1302              : !> \par History
    1303              : !>      05.2005 created [MI]
    1304              : ! **************************************************************************************************
    1305          788 :    SUBROUTINE read_loc_section(localized_wfn_control, loc_section, &
    1306              :                                localize, do_mixed, do_xas, nloc_xas, spin_channel_xas)
    1307              : 
    1308              :       TYPE(localized_wfn_control_type), POINTER          :: localized_wfn_control
    1309              :       TYPE(section_vals_type), POINTER                   :: loc_section
    1310              :       LOGICAL, INTENT(OUT)                               :: localize
    1311              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_mixed, do_xas
    1312              :       INTEGER, INTENT(IN), OPTIONAL                      :: nloc_xas, spin_channel_xas
    1313              : 
    1314              :       INTEGER                                            :: i, ind, ir, n_list, n_rep, n_state, &
    1315              :                                                             nextra, nline, other_spin, &
    1316              :                                                             output_unit, spin_xas
    1317          394 :       INTEGER, DIMENSION(:), POINTER                     :: list, loc_list
    1318              :       LOGICAL                                            :: my_do_mixed, my_do_xas
    1319          394 :       REAL(dp), POINTER                                  :: ene(:)
    1320              :       TYPE(cp_logger_type), POINTER                      :: logger
    1321              :       TYPE(section_vals_type), POINTER                   :: loc_print_section
    1322              : 
    1323          394 :       my_do_xas = .FALSE.
    1324          394 :       spin_xas = 1
    1325          394 :       IF (PRESENT(do_xas)) THEN
    1326          100 :          my_do_xas = do_xas
    1327          100 :          CPASSERT(PRESENT(nloc_xas))
    1328              :       END IF
    1329          394 :       IF (PRESENT(spin_channel_xas)) spin_xas = spin_channel_xas
    1330          394 :       my_do_mixed = .FALSE.
    1331          394 :       IF (PRESENT(do_mixed)) THEN
    1332          394 :          my_do_mixed = do_mixed
    1333              :       END IF
    1334          394 :       CPASSERT(ASSOCIATED(loc_section))
    1335          394 :       NULLIFY (logger)
    1336          394 :       logger => cp_get_default_logger()
    1337              : 
    1338          394 :       CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=localize)
    1339          394 :       IF (localize) THEN
    1340          354 :          loc_print_section => section_vals_get_subs_vals(loc_section, "PRINT")
    1341          354 :          NULLIFY (list)
    1342          354 :          NULLIFY (loc_list)
    1343         2478 :          localized_wfn_control%lu_bound_states = 0
    1344         1062 :          localized_wfn_control%lu_ene_bound = 0.0_dp
    1345         1062 :          localized_wfn_control%nloc_states = 0
    1346          354 :          localized_wfn_control%set_of_states = 0
    1347          354 :          localized_wfn_control%nextra = 0
    1348          354 :          n_state = 0
    1349              : 
    1350              :          CALL section_vals_val_get(loc_section, "MAX_ITER", &
    1351          354 :                                    i_val=localized_wfn_control%max_iter)
    1352              :          CALL section_vals_val_get(loc_section, "MAX_CRAZY_ANGLE", &
    1353          354 :                                    r_val=localized_wfn_control%max_crazy_angle)
    1354              :          CALL section_vals_val_get(loc_section, "CRAZY_SCALE", &
    1355          354 :                                    r_val=localized_wfn_control%crazy_scale)
    1356              :          CALL section_vals_val_get(loc_section, "EPS_OCCUPATION", &
    1357          354 :                                    r_val=localized_wfn_control%eps_occ)
    1358              :          CALL section_vals_val_get(loc_section, "CRAZY_USE_DIAG", &
    1359          354 :                                    l_val=localized_wfn_control%crazy_use_diag)
    1360              :          CALL section_vals_val_get(loc_section, "OUT_ITER_EACH", &
    1361          354 :                                    i_val=localized_wfn_control%out_each)
    1362              :          CALL section_vals_val_get(loc_section, "EPS_LOCALIZATION", &
    1363          354 :                                    r_val=localized_wfn_control%eps_localization)
    1364              :          CALL section_vals_val_get(loc_section, "MIN_OR_MAX", &
    1365          354 :                                    i_val=localized_wfn_control%min_or_max)
    1366              :          CALL section_vals_val_get(loc_section, "JACOBI_FALLBACK", &
    1367          354 :                                    l_val=localized_wfn_control%jacobi_fallback)
    1368              :          CALL section_vals_val_get(loc_section, "JACOBI_REFINEMENT", &
    1369          354 :                                    l_val=localized_wfn_control%jacobi_refinement)
    1370              :          CALL section_vals_val_get(loc_section, "METHOD", &
    1371          354 :                                    i_val=localized_wfn_control%localization_method)
    1372              :          CALL section_vals_val_get(loc_section, "OPERATOR", &
    1373          354 :                                    i_val=localized_wfn_control%operator_type)
    1374              :          CALL section_vals_val_get(loc_section, "RESTART", &
    1375          354 :                                    l_val=localized_wfn_control%loc_restart)
    1376              :          CALL section_vals_val_get(loc_section, "USE_HISTORY", &
    1377          354 :                                    l_val=localized_wfn_control%use_history)
    1378              :          CALL section_vals_val_get(loc_section, "NEXTRA", &
    1379          354 :                                    i_val=localized_wfn_control%nextra)
    1380              :          CALL section_vals_val_get(loc_section, "CPO_GUESS", &
    1381          354 :                                    i_val=localized_wfn_control%coeff_po_guess)
    1382              :          CALL section_vals_val_get(loc_section, "CPO_GUESS_SPACE", &
    1383          354 :                                    i_val=localized_wfn_control%coeff_po_guess_mo_space)
    1384              :          CALL section_vals_val_get(loc_section, "CG_PO", &
    1385          354 :                                    l_val=localized_wfn_control%do_cg_po)
    1386              : 
    1387          354 :          IF (localized_wfn_control%do_homo) THEN
    1388              :             ! List of States HOMO
    1389          346 :             CALL section_vals_val_get(loc_section, "LIST", n_rep_val=n_rep)
    1390          346 :             IF (n_rep > 0) THEN
    1391           14 :                n_list = 0
    1392           28 :                DO ir = 1, n_rep
    1393           14 :                   NULLIFY (list)
    1394           14 :                   CALL section_vals_val_get(loc_section, "LIST", i_rep_val=ir, i_vals=list)
    1395           28 :                   IF (ASSOCIATED(list)) THEN
    1396           14 :                      CALL reallocate(loc_list, 1, n_list + SIZE(list))
    1397           90 :                      DO i = 1, SIZE(list)
    1398           90 :                         loc_list(n_list + i) = list(i)
    1399              :                      END DO ! i
    1400           14 :                      n_list = n_list + SIZE(list)
    1401              :                   END IF
    1402              :                END DO ! ir
    1403           14 :                IF (n_list /= 0) THEN
    1404           14 :                   localized_wfn_control%set_of_states = state_loc_list
    1405           42 :                   ALLOCATE (localized_wfn_control%loc_states(n_list, 2))
    1406          194 :                   localized_wfn_control%loc_states = 0
    1407          180 :                   localized_wfn_control%loc_states(:, 1) = loc_list(:)
    1408          180 :                   localized_wfn_control%loc_states(:, 2) = loc_list(:)
    1409           14 :                   localized_wfn_control%nloc_states(1) = n_list
    1410           14 :                   localized_wfn_control%nloc_states(2) = n_list
    1411           14 :                   IF (my_do_xas) THEN
    1412            4 :                      other_spin = 2
    1413            4 :                      IF (spin_xas == 2) other_spin = 1
    1414            4 :                      localized_wfn_control%nloc_states(other_spin) = 0
    1415           22 :                      localized_wfn_control%loc_states(:, other_spin) = 0
    1416              :                   END IF
    1417           14 :                   DEALLOCATE (loc_list)
    1418              :                END IF
    1419              :             END IF
    1420              : 
    1421              :          ELSE
    1422              :             ! List of States LUMO
    1423            8 :             CALL section_vals_val_get(loc_section, "LIST_UNOCCUPIED", n_rep_val=n_rep)
    1424            8 :             IF (n_rep > 0) THEN
    1425            6 :                n_list = 0
    1426           12 :                DO ir = 1, n_rep
    1427            6 :                   NULLIFY (list)
    1428            6 :                   CALL section_vals_val_get(loc_section, "LIST_UNOCCUPIED", i_rep_val=ir, i_vals=list)
    1429           12 :                   IF (ASSOCIATED(list)) THEN
    1430            6 :                      CALL reallocate(loc_list, 1, n_list + SIZE(list))
    1431           46 :                      DO i = 1, SIZE(list)
    1432           46 :                         loc_list(n_list + i) = list(i)
    1433              :                      END DO ! i
    1434            6 :                      n_list = n_list + SIZE(list)
    1435              :                   END IF
    1436              :                END DO ! ir
    1437            6 :                IF (n_list /= 0) THEN
    1438            6 :                   localized_wfn_control%set_of_states = state_loc_list
    1439           18 :                   ALLOCATE (localized_wfn_control%loc_states(n_list, 2))
    1440           98 :                   localized_wfn_control%loc_states = 0
    1441           92 :                   localized_wfn_control%loc_states(:, 1) = loc_list(:)
    1442           92 :                   localized_wfn_control%loc_states(:, 2) = loc_list(:)
    1443            6 :                   localized_wfn_control%nloc_states(1) = n_list
    1444            6 :                   DEALLOCATE (loc_list)
    1445              :                END IF
    1446              :             END IF
    1447              :          END IF
    1448              : 
    1449          354 :          IF (localized_wfn_control%set_of_states == 0) THEN
    1450          334 :             CALL section_vals_val_get(loc_section, "ENERGY_RANGE", r_vals=ene)
    1451          334 :             IF (ene(1) /= ene(2)) THEN
    1452           10 :                localized_wfn_control%set_of_states = energy_loc_range
    1453           10 :                localized_wfn_control%lu_ene_bound(1) = ene(1)
    1454           10 :                localized_wfn_control%lu_ene_bound(2) = ene(2)
    1455              :             END IF
    1456              :          END IF
    1457              : 
    1458              :          ! All States or XAS specific states
    1459          354 :          IF (localized_wfn_control%set_of_states == 0) THEN
    1460          324 :             IF (my_do_xas) THEN
    1461           72 :                localized_wfn_control%set_of_states = state_loc_range
    1462          216 :                localized_wfn_control%nloc_states(:) = 0
    1463          216 :                localized_wfn_control%lu_bound_states(1, :) = 0
    1464          216 :                localized_wfn_control%lu_bound_states(2, :) = 0
    1465           72 :                localized_wfn_control%nloc_states(spin_xas) = nloc_xas
    1466           72 :                localized_wfn_control%lu_bound_states(1, spin_xas) = 1
    1467           72 :                localized_wfn_control%lu_bound_states(2, spin_xas) = nloc_xas
    1468          252 :             ELSE IF (my_do_mixed) THEN
    1469            2 :                localized_wfn_control%set_of_states = state_loc_mixed
    1470            2 :                nextra = localized_wfn_control%nextra
    1471              :             ELSE
    1472          250 :                localized_wfn_control%set_of_states = state_loc_all
    1473              :             END IF
    1474              :          END IF
    1475              : 
    1476              :          localized_wfn_control%print_centers = &
    1477              :             BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
    1478          354 :                                              "WANNIER_CENTERS"), cp_p_file)
    1479              :          localized_wfn_control%print_spreads = &
    1480              :             BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
    1481          354 :                                              "WANNIER_SPREADS"), cp_p_file)
    1482              :          localized_wfn_control%print_cubes = &
    1483              :             BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
    1484          354 :                                              "WANNIER_CUBES"), cp_p_file)
    1485              : 
    1486              :          output_unit = cp_print_key_unit_nr(logger, loc_print_section, "PROGRAM_RUN_INFO", &
    1487          354 :                                             extension=".Log")
    1488              : 
    1489          354 :          IF (output_unit > 0) THEN
    1490              :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
    1491          177 :                "LOCALIZE| The spread relative to a set of orbitals is computed"
    1492              : 
    1493          302 :             SELECT CASE (localized_wfn_control%set_of_states)
    1494              :             CASE (state_loc_all)
    1495              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1496          125 :                   "LOCALIZE| Orbitals to be localized: All orbitals"
    1497              :                WRITE (UNIT=output_unit, FMT="(T2,A,/,T12,A,F16.8)") &
    1498          125 :                   "LOCALIZE| If fractional occupation, fully occupied MOs are those ", &
    1499          250 :                   "within occupation tolerance of ", localized_wfn_control%eps_occ
    1500              :             CASE (state_loc_range)
    1501              :                WRITE (UNIT=output_unit, FMT="(T2,A,T65,I8,A,I8)") &
    1502           36 :                   "LOCALIZE| Orbitals to be localized: Those with index between ", &
    1503           36 :                   localized_wfn_control%lu_bound_states(1, spin_xas), " and ", &
    1504           72 :                   localized_wfn_control%lu_bound_states(2, spin_xas)
    1505              :             CASE (state_loc_list)
    1506              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1507           10 :                   "LOCALIZE| Orbitals to be localized: Those with index in the following list"
    1508           10 :                nline = localized_wfn_control%nloc_states(1)/10 + 1
    1509           10 :                ind = 0
    1510           21 :                DO i = 1, nline
    1511           21 :                   IF (ind + 10 < localized_wfn_control%nloc_states(1)) THEN
    1512            1 :                      WRITE (UNIT=output_unit, FMT="(T8,10I7)") localized_wfn_control%loc_states(ind + 1:ind + 10, 1)
    1513            1 :                      ind = ind + 10
    1514              :                   ELSE
    1515              :                      WRITE (UNIT=output_unit, FMT="(T8,10I7)") &
    1516           10 :                         localized_wfn_control%loc_states(ind + 1:localized_wfn_control%nloc_states(1), 1)
    1517           10 :                      ind = localized_wfn_control%nloc_states(1)
    1518              :                   END IF
    1519              :                END DO
    1520              :             CASE (energy_loc_range)
    1521              :                WRITE (UNIT=output_unit, FMT="(T2,A,T65,/,f16.6,A,f16.6,A)") &
    1522            5 :                   "LOCALIZE| Orbitals to be localized: Those with energy in the range between ", &
    1523           10 :                   localized_wfn_control%lu_ene_bound(1), " and ", localized_wfn_control%lu_ene_bound(2), " a.u."
    1524              :             CASE (state_loc_mixed)
    1525              :                WRITE (UNIT=output_unit, FMT="(T2,A,I4,A)") &
    1526            1 :                   "LOCALIZE| Orbitals to be localized: Occupied orbitals + ", nextra, " orbitals"
    1527              :             CASE DEFAULT
    1528              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1529          177 :                   "LOCALIZE| Orbitals to be localized: None "
    1530              :             END SELECT
    1531              : 
    1532          354 :             SELECT CASE (localized_wfn_control%operator_type)
    1533              :             CASE (op_loc_berry)
    1534              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1535          177 :                   "LOCALIZE| Spread defined by the Berry phase operator "
    1536              :             CASE (op_loc_boys)
    1537              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1538            0 :                   "LOCALIZE| Spread defined by the Boys phase operator "
    1539              :             CASE DEFAULT
    1540              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1541          177 :                   "LOCALIZE| Spread defined by the Pipek phase operator "
    1542              :             END SELECT
    1543              : 
    1544          313 :             SELECT CASE (localized_wfn_control%localization_method)
    1545              :             CASE (do_loc_jacobi)
    1546              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1547          136 :                   "LOCALIZE| Optimal unitary transformation generated by Jacobi algorithm"
    1548              :             CASE (do_loc_crazy)
    1549              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1550           29 :                   "LOCALIZE| Optimal unitary transformation generated by Crazy angle algorithm"
    1551              :                WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
    1552           29 :                   "LOCALIZE| maximum angle: ", localized_wfn_control%max_crazy_angle
    1553              :                WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
    1554           29 :                   "LOCALIZE| scaling: ", localized_wfn_control%crazy_scale
    1555              :                WRITE (UNIT=output_unit, FMT="(T2,A,L1)") &
    1556           29 :                   "LOCALIZE| use diag:", localized_wfn_control%crazy_use_diag
    1557              :             CASE (do_loc_gapo)
    1558              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1559            1 :                   "LOCALIZE| Optimal unitary transformation generated by gradient ascent algorithm "
    1560              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1561            1 :                   "LOCALIZE| for partially occupied wannier functions"
    1562              :             CASE (do_loc_direct)
    1563              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1564            1 :                   "LOCALIZE| Optimal unitary transformation generated by direct algorithm"
    1565              :             CASE (do_loc_l1_norm_sd)
    1566              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1567            9 :                   "LOCALIZE| Optimal unitary transformation generated by "
    1568              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1569            9 :                   "LOCALIZE| steepest descent algorithm applied on an approximate l1 norm"
    1570              :             CASE (do_loc_none)
    1571              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1572            0 :                   "LOCALIZE| No unitary transformation is applied"
    1573              :             CASE (do_loc_scdm)
    1574              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    1575          177 :                   "LOCALIZE| Pivoted QR decomposition is used to transform coefficients"
    1576              :             END SELECT
    1577              : 
    1578              :          END IF ! process has output_unit
    1579              : 
    1580          354 :          CALL cp_print_key_finished_output(output_unit, logger, loc_print_section, "PROGRAM_RUN_INFO")
    1581              : 
    1582              :       ELSE
    1583           40 :          localized_wfn_control%localization_method = do_loc_none
    1584           40 :          localized_wfn_control%localization_method = state_loc_none
    1585           40 :          localized_wfn_control%print_centers = .FALSE.
    1586           40 :          localized_wfn_control%print_spreads = .FALSE.
    1587           40 :          localized_wfn_control%print_cubes = .FALSE.
    1588              :       END IF
    1589              : 
    1590          394 :    END SUBROUTINE read_loc_section
    1591              : 
    1592              : ! **************************************************************************************************
    1593              : !> \brief create the center and spread array and the file names for the output
    1594              : !> \param localized_wfn_control ...
    1595              : !> \param nmoloc ...
    1596              : !> \param nspins ...
    1597              : !> \par History
    1598              : !>      04.2005 created [MI]
    1599              : ! **************************************************************************************************
    1600          446 :    SUBROUTINE set_loc_centers(localized_wfn_control, nmoloc, nspins)
    1601              : 
    1602              :       TYPE(localized_wfn_control_type)                   :: localized_wfn_control
    1603              :       INTEGER, DIMENSION(2), INTENT(IN)                  :: nmoloc
    1604              :       INTEGER, INTENT(IN)                                :: nspins
    1605              : 
    1606              :       INTEGER                                            :: ispin
    1607              : 
    1608         1046 :       DO ispin = 1, nspins
    1609         1758 :          ALLOCATE (localized_wfn_control%centers_set(ispin)%array(6, nmoloc(ispin)))
    1610        28150 :          localized_wfn_control%centers_set(ispin)%array = 0.0_dp
    1611              :       END DO
    1612              : 
    1613          446 :    END SUBROUTINE set_loc_centers
    1614              : 
    1615              : ! **************************************************************************************************
    1616              : !> \brief create the lists of mos that are taken into account
    1617              : !> \param localized_wfn_control ...
    1618              : !> \param nmoloc ...
    1619              : !> \param nmo ...
    1620              : !> \param nspins ...
    1621              : !> \param my_spin ...
    1622              : !> \par History
    1623              : !>      04.2005 created [MI]
    1624              : ! **************************************************************************************************
    1625          316 :    SUBROUTINE set_loc_wfn_lists(localized_wfn_control, nmoloc, nmo, nspins, my_spin)
    1626              : 
    1627              :       TYPE(localized_wfn_control_type)                   :: localized_wfn_control
    1628              :       INTEGER, DIMENSION(2), INTENT(IN)                  :: nmoloc, nmo
    1629              :       INTEGER, INTENT(IN)                                :: nspins
    1630              :       INTEGER, INTENT(IN), OPTIONAL                      :: my_spin
    1631              : 
    1632              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'set_loc_wfn_lists'
    1633              : 
    1634              :       INTEGER                                            :: i, ispin, max_iloc, max_nmoloc, state
    1635              : 
    1636          316 :       CALL timeset(routineN, state)
    1637              : 
    1638          948 :       localized_wfn_control%nloc_states(1:2) = nmoloc(1:2)
    1639          316 :       max_nmoloc = MAX(nmoloc(1), nmoloc(2))
    1640              : 
    1641          334 :       SELECT CASE (localized_wfn_control%set_of_states)
    1642              :       CASE (state_loc_list)
    1643              :          ! List
    1644           18 :          CPASSERT(ASSOCIATED(localized_wfn_control%loc_states))
    1645           52 :          DO ispin = 1, nspins
    1646           34 :             localized_wfn_control%lu_bound_states(1, ispin) = 1
    1647           34 :             localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
    1648           52 :             IF (nmoloc(ispin) < 1) THEN
    1649            4 :                localized_wfn_control%lu_bound_states(1, ispin) = 0
    1650           22 :                localized_wfn_control%loc_states(:, ispin) = 0
    1651              :             END IF
    1652              :          END DO
    1653              :       CASE (state_loc_range)
    1654              :          ! Range
    1655          114 :          ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
    1656          446 :          localized_wfn_control%loc_states = 0
    1657          114 :          DO ispin = 1, nspins
    1658              :             localized_wfn_control%lu_bound_states(1, ispin) = &
    1659           76 :                localized_wfn_control%lu_bound_states(1, my_spin)
    1660              :             localized_wfn_control%lu_bound_states(2, ispin) = &
    1661           76 :                localized_wfn_control%lu_bound_states(1, my_spin) + nmoloc(ispin) - 1
    1662           76 :             max_iloc = localized_wfn_control%lu_bound_states(2, ispin)
    1663          242 :             DO i = 1, nmoloc(ispin)
    1664          242 :                localized_wfn_control%loc_states(i, ispin) = localized_wfn_control%lu_bound_states(1, ispin) + i - 1
    1665              :             END DO
    1666           76 :             CPASSERT(max_iloc <= nmo(ispin))
    1667           38 :             MARK_USED(nmo)
    1668              :          END DO
    1669              :       CASE (energy_loc_range)
    1670              :          ! Energy
    1671           30 :          ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
    1672          214 :          localized_wfn_control%loc_states = 0
    1673           22 :          DO ispin = 1, nspins
    1674          134 :             DO i = 1, nmoloc(ispin)
    1675          124 :                localized_wfn_control%loc_states(i, ispin) = localized_wfn_control%lu_bound_states(1, ispin) + i - 1
    1676              :             END DO
    1677              :          END DO
    1678              :       CASE (state_loc_all)
    1679              :          ! All
    1680          744 :          ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
    1681         4624 :          localized_wfn_control%loc_states = 0
    1682              : 
    1683          248 :          IF (localized_wfn_control%lu_bound_states(1, 1) == 1) THEN
    1684          582 :             DO ispin = 1, nspins
    1685          334 :                localized_wfn_control%lu_bound_states(1, ispin) = 1
    1686          334 :                localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
    1687          334 :                IF (nmoloc(ispin) < 1) localized_wfn_control%lu_bound_states(1, ispin) = 0
    1688         3176 :                DO i = 1, nmoloc(ispin)
    1689         2928 :                   localized_wfn_control%loc_states(i, ispin) = i
    1690              :                END DO
    1691              :             END DO
    1692              :          ELSE
    1693            0 :             DO ispin = 1, nspins
    1694            0 :                IF (nmoloc(ispin) < 1) localized_wfn_control%lu_bound_states(1, ispin) = 0
    1695            0 :                DO i = 1, nmoloc(ispin)
    1696              :                   localized_wfn_control%loc_states(i, ispin) = &
    1697            0 :                      localized_wfn_control%lu_bound_states(1, ispin) + i - 1
    1698              :                END DO
    1699              :             END DO
    1700              :          END IF
    1701              :       CASE (state_loc_mixed)
    1702              :          ! Mixed
    1703            6 :          ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
    1704          178 :          localized_wfn_control%loc_states = 0
    1705          320 :          DO ispin = 1, nspins
    1706           90 :             DO i = 1, nmoloc(ispin)
    1707           88 :                localized_wfn_control%loc_states(i, ispin) = i
    1708              :             END DO
    1709              :          END DO
    1710              :       END SELECT
    1711              : 
    1712          316 :       CALL timestop(state)
    1713              : 
    1714          316 :    END SUBROUTINE set_loc_wfn_lists
    1715              : 
    1716              : END MODULE qs_loc_utils
    1717              : 
        

Generated by: LCOV version 2.0-1