LCOV - code coverage report
Current view: top level - src - almo_scf_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 813 1068 76.1 %
Date: 2024-03-28 07:31:50 Functions: 20 25 80.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Subroutines for ALMO SCF
      10             : !> \par History
      11             : !>       2011.06 created [Rustam Z Khaliullin]
      12             : !>       2018.09 smearing support [Ruben Staub]
      13             : !> \author Rustam Z Khaliullin
      14             : ! **************************************************************************************************
      15             : MODULE almo_scf_methods
      16             :    USE almo_scf_types,                  ONLY: almo_scf_env_type,&
      17             :                                               almo_scf_history_type
      18             :    USE bibliography,                    ONLY: Kolafa2004,&
      19             :                                               Kuhne2007,&
      20             :                                               cite_reference
      21             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      22             :    USE cp_dbcsr_cholesky,               ONLY: cp_dbcsr_cholesky_decompose,&
      23             :                                               cp_dbcsr_cholesky_invert
      24             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      25             :                                               cp_logger_get_default_unit_nr,&
      26             :                                               cp_logger_type
      27             :    USE dbcsr_api,                       ONLY: &
      28             :         dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, &
      29             :         dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_frobenius_norm, &
      30             :         dbcsr_get_block_p, dbcsr_get_diag, dbcsr_get_info, dbcsr_get_stored_coordinates, &
      31             :         dbcsr_init_random, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      32             :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
      33             :         dbcsr_nblkcols_total, dbcsr_nblkrows_total, dbcsr_print, dbcsr_release, &
      34             :         dbcsr_reserve_block2d, dbcsr_scale_by_vector, dbcsr_set, dbcsr_set_diag, dbcsr_transposed, &
      35             :         dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric, dbcsr_work_create
      36             :    USE domain_submatrix_methods,        ONLY: &
      37             :         add_submatrices, construct_dbcsr_from_submatrices, construct_submatrices, &
      38             :         copy_submatrices, copy_submatrix_data, init_submatrices, multiply_submatrices, &
      39             :         print_submatrices, release_submatrices
      40             :    USE domain_submatrix_types,          ONLY: domain_map_type,&
      41             :                                               domain_submatrix_type,&
      42             :                                               select_row,&
      43             :                                               select_row_col
      44             :    USE fermi_utils,                     ONLY: FermiFixed
      45             : !! for smearing
      46             :    USE input_constants,                 ONLY: almo_domain_layout_molecular,&
      47             :                                               almo_mat_distr_atomic,&
      48             :                                               almo_scf_diag,&
      49             :                                               spd_inversion_dense_cholesky,&
      50             :                                               spd_inversion_ls_hotelling,&
      51             :                                               spd_inversion_ls_taylor
      52             :    USE iterate_matrix,                  ONLY: invert_Hotelling,&
      53             :                                               invert_Taylor,&
      54             :                                               matrix_sqrt_Newton_Schulz
      55             :    USE kinds,                           ONLY: dp
      56             :    USE mathlib,                         ONLY: binomial
      57             :    USE message_passing,                 ONLY: mp_comm_type,&
      58             :                                               mp_para_env_type
      59             :    USE util,                            ONLY: sort
      60             : #include "./base/base_uses.f90"
      61             : 
      62             :    IMPLICIT NONE
      63             : 
      64             :    PRIVATE
      65             : 
      66             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf_methods'
      67             : 
      68             :    PUBLIC almo_scf_ks_to_ks_blk, almo_scf_p_blk_to_t_blk, &
      69             :       almo_scf_t_to_proj, almo_scf_ks_blk_to_tv_blk, &
      70             :       almo_scf_ks_xx_to_tv_xx, almo_scf_t_rescaling, &
      71             :       apply_projector, get_overlap, &
      72             :       generator_to_unitary, &
      73             :       orthogonalize_mos, &
      74             :       pseudo_invert_diagonal_blk, construct_test, &
      75             :       construct_domain_preconditioner, &
      76             :       apply_domain_operators, &
      77             :       construct_domain_s_inv, &
      78             :       construct_domain_s_sqrt, &
      79             :       distribute_domains, &
      80             :       almo_scf_ks_to_ks_xx, &
      81             :       construct_domain_r_down, &
      82             :       xalmo_initial_guess, &
      83             :       fill_matrix_with_ones
      84             : 
      85             : CONTAINS
      86             : 
      87             : ! **************************************************************************************************
      88             : !> \brief Fill all matrix blocks with 1.0_dp
      89             : !> \param matrix ...
      90             : !> \par History
      91             : !>       2019.09 created [Rustam Z Khaliullin]
      92             : !> \author Rustam Z Khaliullin
      93             : ! **************************************************************************************************
      94           0 :    SUBROUTINE fill_matrix_with_ones(matrix)
      95             : 
      96             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
      97             : 
      98             :       INTEGER                                            :: col, hold, iblock_col, iblock_row, &
      99             :                                                             mynode, nblkcols_tot, nblkrows_tot, row
     100             :       LOGICAL                                            :: tr
     101           0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_new_block
     102             :       TYPE(dbcsr_distribution_type)                      :: dist
     103             : 
     104           0 :       CALL dbcsr_get_info(matrix, distribution=dist)
     105           0 :       CALL dbcsr_distribution_get(dist, mynode=mynode)
     106           0 :       CALL dbcsr_work_create(matrix, work_mutable=.TRUE.)
     107           0 :       nblkrows_tot = dbcsr_nblkrows_total(matrix)
     108           0 :       nblkcols_tot = dbcsr_nblkcols_total(matrix)
     109           0 :       DO row = 1, nblkrows_tot
     110           0 :          DO col = 1, nblkcols_tot
     111           0 :             tr = .FALSE.
     112           0 :             iblock_row = row
     113           0 :             iblock_col = col
     114             :             CALL dbcsr_get_stored_coordinates(matrix, &
     115           0 :                                               iblock_row, iblock_col, hold)
     116           0 :             IF (hold .EQ. mynode) THEN
     117           0 :                NULLIFY (p_new_block)
     118             :                CALL dbcsr_reserve_block2d(matrix, &
     119           0 :                                           iblock_row, iblock_col, p_new_block)
     120           0 :                CPASSERT(ASSOCIATED(p_new_block))
     121           0 :                p_new_block(:, :) = 1.0_dp
     122             :             END IF
     123             :          END DO
     124             :       END DO
     125           0 :       CALL dbcsr_finalize(matrix)
     126             : 
     127           0 :    END SUBROUTINE fill_matrix_with_ones
     128             : 
     129             : ! **************************************************************************************************
     130             : !> \brief builds projected KS matrices for the overlapping domains
     131             : !>        also computes the DIIS error vector as a by-product
     132             : !> \param almo_scf_env ...
     133             : !> \par History
     134             : !>       2013.03 created [Rustam Z Khaliullin]
     135             : !> \author Rustam Z Khaliullin
     136             : ! **************************************************************************************************
     137           2 :    SUBROUTINE almo_scf_ks_to_ks_xx(almo_scf_env)
     138             : 
     139             :       TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
     140             : 
     141             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'almo_scf_ks_to_ks_xx'
     142             : 
     143             :       INTEGER                                            :: handle, ispin, ndomains
     144             :       REAL(KIND=dp)                                      :: eps_multiply
     145             :       TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2, matrix_tmp3, matrix_tmp4, matrix_tmp5, &
     146             :          matrix_tmp6, matrix_tmp7, matrix_tmp8, matrix_tmp9
     147             :       TYPE(domain_submatrix_type), ALLOCATABLE, &
     148           2 :          DIMENSION(:)                                    :: subm_tmp1, subm_tmp2, subm_tmp3
     149             : 
     150           2 :       CALL timeset(routineN, handle)
     151             : 
     152           2 :       eps_multiply = almo_scf_env%eps_filter
     153             : 
     154           4 :       DO ispin = 1, almo_scf_env%nspins
     155             : 
     156           2 :          ndomains = dbcsr_nblkcols_total(almo_scf_env%quench_t(ispin))
     157             : 
     158             :          ! 0. Create KS_xx
     159             :          CALL construct_submatrices( &
     160             :             almo_scf_env%matrix_ks(ispin), &
     161             :             almo_scf_env%domain_ks_xx(:, ispin), &
     162             :             almo_scf_env%quench_t(ispin), &
     163             :             almo_scf_env%domain_map(ispin), &
     164             :             almo_scf_env%cpu_of_domain, &
     165           2 :             select_row_col)
     166             : 
     167             :          !!!!! RZK-warning MAKE SURE THAT YOU NEED BLOCKS OUTSIDE QUENCH_T
     168             :          !!!!! FOR ALL NO-MATRICES NOT COMPUTING THEM CAN SAVE LOTS OF TIME
     169             : 
     170             :          ! 1. TMP1=KS.T
     171             :          !    Cost: NOn
     172             :          !matrix_tmp1 = create NxO, full
     173             :          CALL dbcsr_create(matrix_tmp1, &
     174           2 :                            template=almo_scf_env%matrix_t(ispin))
     175             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_ks(ispin), &
     176             :                              almo_scf_env%matrix_t(ispin), &
     177             :                              0.0_dp, matrix_tmp1, &
     178           2 :                              filter_eps=eps_multiply)
     179             : 
     180             :          ! 2. TMP2=TMP1.SigInv=KS.T.SigInv
     181             :          !    Cost: NOO
     182             :          !matrix_tmp2 = create NxO, full
     183             :          CALL dbcsr_create(matrix_tmp2, &
     184           2 :                            template=almo_scf_env%matrix_t(ispin))
     185             :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, &
     186             :                              almo_scf_env%matrix_sigma_inv(ispin), &
     187             :                              0.0_dp, matrix_tmp2, &
     188           2 :                              filter_eps=eps_multiply)
     189             : 
     190             :          ! 3. TMP1=S.T
     191             :          !    Cost: NOn
     192             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
     193             :                              almo_scf_env%matrix_t(ispin), &
     194             :                              0.0_dp, matrix_tmp1, &
     195           2 :                              filter_eps=eps_multiply)
     196             : 
     197             :          ! 4. TMP4=TMP2.tr(TMP1)=KS.T.SigInv.tr(T).S
     198             :          !    Cost: NNO
     199             :          !matrix_tmp4 = create NxN
     200             :          CALL dbcsr_create(matrix_tmp4, &
     201             :                            template=almo_scf_env%matrix_s(1), &
     202           2 :                            matrix_type=dbcsr_type_no_symmetry)
     203             :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp2, &
     204             :                              matrix_tmp1, &
     205             :                              0.0_dp, matrix_tmp4, &
     206           2 :                              filter_eps=eps_multiply)
     207             : 
     208             :          ! 5. KS_xx=KS_xx-TMP4_xx-tr(TMP4_xx)
     209          16 :          ALLOCATE (subm_tmp1(ndomains))
     210           2 :          CALL init_submatrices(subm_tmp1)
     211             :          CALL construct_submatrices( &
     212             :             matrix_tmp4, &
     213             :             subm_tmp1, &
     214             :             almo_scf_env%quench_t(ispin), &
     215             :             almo_scf_env%domain_map(ispin), &
     216             :             almo_scf_env%cpu_of_domain, &
     217           2 :             select_row_col)
     218             :          CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
     219           2 :                               -1.0_dp, subm_tmp1, 'N')
     220             :          CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
     221           2 :                               -1.0_dp, subm_tmp1, 'T')
     222             : 
     223             :          ! 6. TMP3=tr(TMP4).T=S.T.SigInv.tr(T).KS.T
     224             :          !    Cost: NOn
     225             :          !matrix_tmp3 = create NxO, full
     226             :          CALL dbcsr_create(matrix_tmp3, &
     227             :                            template=almo_scf_env%matrix_t(ispin), &
     228           2 :                            matrix_type=dbcsr_type_no_symmetry)
     229             :          CALL dbcsr_multiply("T", "N", 1.0_dp, &
     230             :                              matrix_tmp4, &
     231             :                              almo_scf_env%matrix_t(ispin), &
     232             :                              0.0_dp, matrix_tmp3, &
     233           2 :                              filter_eps=eps_multiply)
     234           2 :          CALL dbcsr_release(matrix_tmp4)
     235             : 
     236             :          ! 8. TMP6=TMP3.SigInv=S.T.SigInv.tr(T).KS.T.SigInv
     237             :          !    Cost: NOO
     238             :          !matrix_tmp6 = create NxO, full
     239             :          CALL dbcsr_create(matrix_tmp6, &
     240             :                            template=almo_scf_env%matrix_t(ispin), &
     241           2 :                            matrix_type=dbcsr_type_no_symmetry)
     242             :          CALL dbcsr_multiply("N", "N", 1.0_dp, &
     243             :                              matrix_tmp3, &
     244             :                              almo_scf_env%matrix_sigma_inv(ispin), &
     245             :                              0.0_dp, matrix_tmp6, &
     246           2 :                              filter_eps=eps_multiply)
     247             : 
     248             :          ! 8A. Use intermediate matrices to evaluate the gradient/error
     249             :          !     Err=(TMP2-TMP6)_q=(KS.T.SigInv-S.T.SigInv.tr(T).KS.T.SigInv)_q
     250             :          ! error vector in AO-MO basis
     251             :          CALL dbcsr_copy(almo_scf_env%matrix_err_xx(ispin), &
     252           2 :                          almo_scf_env%quench_t(ispin))
     253             :          CALL dbcsr_copy(almo_scf_env%matrix_err_xx(ispin), &
     254           2 :                          matrix_tmp2, keep_sparsity=.TRUE.)
     255             :          CALL dbcsr_create(matrix_tmp4, &
     256             :                            template=almo_scf_env%matrix_t(ispin), &
     257           2 :                            matrix_type=dbcsr_type_no_symmetry)
     258             :          CALL dbcsr_copy(matrix_tmp4, &
     259           2 :                          almo_scf_env%quench_t(ispin))
     260             :          CALL dbcsr_copy(matrix_tmp4, &
     261           2 :                          matrix_tmp6, keep_sparsity=.TRUE.)
     262             :          CALL dbcsr_add(almo_scf_env%matrix_err_xx(ispin), &
     263           2 :                         matrix_tmp4, 1.0_dp, -1.0_dp)
     264           2 :          CALL dbcsr_release(matrix_tmp4)
     265             :          !
     266             :          ! error vector in AO-AO basis
     267             :          ! RZK-warning tmp4 can be created using the sparsity pattern,
     268             :          ! then retain_sparsity can be used to perform the multiply
     269             :          ! this will save some time
     270             :          CALL dbcsr_copy(matrix_tmp3, &
     271           2 :                          matrix_tmp2)
     272             :          CALL dbcsr_add(matrix_tmp3, &
     273           2 :                         matrix_tmp6, 1.0_dp, -1.0_dp)
     274             :          CALL dbcsr_create(matrix_tmp4, &
     275             :                            template=almo_scf_env%matrix_s(1), &
     276           2 :                            matrix_type=dbcsr_type_no_symmetry)
     277             :          CALL dbcsr_multiply("N", "T", 1.0_dp, &
     278             :                              matrix_tmp3, &
     279             :                              almo_scf_env%matrix_t(ispin), &
     280             :                              0.0_dp, matrix_tmp4, &
     281           2 :                              filter_eps=eps_multiply)
     282             :          CALL construct_submatrices( &
     283             :             matrix_tmp4, &
     284             :             almo_scf_env%domain_err(:, ispin), &
     285             :             almo_scf_env%quench_t(ispin), &
     286             :             almo_scf_env%domain_map(ispin), &
     287             :             almo_scf_env%cpu_of_domain, &
     288           2 :             select_row_col)
     289           2 :          CALL dbcsr_release(matrix_tmp4)
     290             :          ! domain_err submatrices are in down-up representation
     291             :          ! bring them into the orthogonalized basis
     292          16 :          ALLOCATE (subm_tmp2(ndomains))
     293           2 :          CALL init_submatrices(subm_tmp2)
     294             :          CALL multiply_submatrices('N', 'N', 1.0_dp, &
     295             :                                    almo_scf_env%domain_err(:, ispin), &
     296           2 :                                    almo_scf_env%domain_s_sqrt(:, ispin), 0.0_dp, subm_tmp2)
     297             :          CALL multiply_submatrices('N', 'N', 1.0_dp, &
     298             :                                    almo_scf_env%domain_s_sqrt_inv(:, ispin), &
     299           2 :                                    subm_tmp2, 0.0_dp, almo_scf_env%domain_err(:, ispin))
     300             : 
     301             :          ! 9. TMP5=TMP6.tr(TMP1)=S.T.SigInv.tr(T).KS.T.SigInv.tr(T).S
     302             :          !    Cost: NNO
     303             :          !    matrix_tmp5 = create NxN, full
     304             :          ! RZK-warning tmp5 can be created using the sparsity pattern,
     305             :          ! then retain_sparsity can be used to perform the multiply
     306             :          ! this will save some time
     307             :          CALL dbcsr_create(matrix_tmp5, &
     308             :                            template=almo_scf_env%matrix_s(1), &
     309           2 :                            matrix_type=dbcsr_type_no_symmetry)
     310             :          CALL dbcsr_multiply("N", "T", 1.0_dp, &
     311             :                              matrix_tmp6, &
     312             :                              matrix_tmp1, &
     313             :                              0.0_dp, matrix_tmp5, &
     314           2 :                              filter_eps=eps_multiply)
     315             : 
     316             :          ! 10. KS_xx=KS_xx+TMP5_xx
     317             :          CALL construct_submatrices( &
     318             :             matrix_tmp5, &
     319             :             subm_tmp1, &
     320             :             almo_scf_env%quench_t(ispin), &
     321             :             almo_scf_env%domain_map(ispin), &
     322             :             almo_scf_env%cpu_of_domain, &
     323           2 :             select_row_col)
     324           2 :          CALL dbcsr_release(matrix_tmp5)
     325             :          CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
     326           2 :                               1.0_dp, subm_tmp1, 'N')
     327             : 
     328             :          ! 11. KS_xx=KS_xx + [S.T]_xx.[SigInv.tr(T).KS.(1-T.SigInv.tr(T).S)]_xx + transposed
     329          16 :          ALLOCATE (subm_tmp3(ndomains))
     330           2 :          CALL init_submatrices(subm_tmp3)
     331             :          CALL construct_submatrices( &
     332             :             matrix_tmp2, &
     333             :             subm_tmp2, &
     334             :             almo_scf_env%quench_t(ispin), &
     335             :             almo_scf_env%domain_map(ispin), &
     336             :             almo_scf_env%cpu_of_domain, &
     337           2 :             select_row)
     338             :          CALL construct_submatrices( &
     339             :             matrix_tmp6, &
     340             :             subm_tmp3, &
     341             :             almo_scf_env%quench_t(ispin), &
     342             :             almo_scf_env%domain_map(ispin), &
     343             :             almo_scf_env%cpu_of_domain, &
     344           2 :             select_row)
     345           2 :          CALL dbcsr_release(matrix_tmp6)
     346             :          CALL add_submatrices(1.0_dp, subm_tmp2, &
     347           2 :                               -1.0_dp, subm_tmp3, 'N')
     348             :          CALL construct_submatrices( &
     349             :             matrix_tmp1, &
     350             :             subm_tmp3, &
     351             :             almo_scf_env%quench_t(ispin), &
     352             :             almo_scf_env%domain_map(ispin), &
     353             :             almo_scf_env%cpu_of_domain, &
     354           2 :             select_row)
     355             :          CALL multiply_submatrices('N', 'T', 1.0_dp, subm_tmp2, &
     356           2 :                                    subm_tmp3, 0.0_dp, subm_tmp1)
     357             :          CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
     358           2 :                               1.0_dp, subm_tmp1, 'N')
     359             :          CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
     360           2 :                               1.0_dp, subm_tmp1, 'T')
     361             : 
     362             :          ! 12. TMP7=tr(T).KS.T.SigInv
     363             :          CALL dbcsr_create(matrix_tmp7, &
     364             :                            template=almo_scf_env%matrix_sigma_blk(ispin), &
     365           2 :                            matrix_type=dbcsr_type_no_symmetry)
     366             :          CALL dbcsr_multiply("T", "N", 1.0_dp, &
     367             :                              almo_scf_env%matrix_t(ispin), &
     368             :                              matrix_tmp2, &
     369             :                              0.0_dp, matrix_tmp7, &
     370           2 :                              filter_eps=eps_multiply)
     371             : 
     372             :          ! 13. TMP8=[SigInv.tr(T).KS.T.SigInv]_xx
     373             :          CALL dbcsr_create(matrix_tmp8, &
     374             :                            template=almo_scf_env%matrix_sigma_blk(ispin), &
     375           2 :                            matrix_type=dbcsr_type_symmetric)
     376           2 :          CALL dbcsr_copy(matrix_tmp8, almo_scf_env%matrix_sigma_blk(ispin))
     377             :          CALL dbcsr_multiply("N", "N", 1.0_dp, &
     378             :                              almo_scf_env%matrix_sigma_inv(ispin), &
     379             :                              matrix_tmp7, &
     380             :                              0.0_dp, matrix_tmp8, &
     381             :                              retain_sparsity=.TRUE., &
     382           2 :                              filter_eps=eps_multiply)
     383           2 :          CALL dbcsr_release(matrix_tmp7)
     384             : 
     385             :          ! 13. TMP9=[S.T]_xx
     386             :          CALL dbcsr_create(matrix_tmp9, &
     387             :                            template=almo_scf_env%matrix_t(ispin), &
     388           2 :                            matrix_type=dbcsr_type_no_symmetry)
     389           2 :          CALL dbcsr_copy(matrix_tmp9, almo_scf_env%quench_t(ispin))
     390           2 :          CALL dbcsr_copy(matrix_tmp9, matrix_tmp1, keep_sparsity=.TRUE.)
     391             : 
     392             :          ! 14. TMP3=TMP9.TMP8=[S.T]_xx.[SigInv.tr(T).KS.T.SigInv]_xx
     393             :          CALL dbcsr_multiply("N", "N", 1.0_dp, &
     394             :                              matrix_tmp9, &
     395             :                              matrix_tmp8, &
     396             :                              0.0_dp, matrix_tmp3, &
     397           2 :                              filter_eps=eps_multiply)
     398           2 :          CALL dbcsr_release(matrix_tmp8)
     399           2 :          CALL dbcsr_release(matrix_tmp9)
     400             : 
     401             :          ! 15. KS_xx=KS_xx+[S.T]_xx.[SigInv.tr(T).KS.T.SigInv]_xx.[tr(T).S]_xx
     402             :          CALL construct_submatrices( &
     403             :             matrix_tmp3, &
     404             :             subm_tmp2, &
     405             :             almo_scf_env%quench_t(ispin), &
     406             :             almo_scf_env%domain_map(ispin), &
     407             :             almo_scf_env%cpu_of_domain, &
     408           2 :             select_row)
     409             :          CALL multiply_submatrices('N', 'T', 1.0_dp, subm_tmp2, &
     410           2 :                                    subm_tmp3, 0.0_dp, subm_tmp1)
     411             :          CALL add_submatrices(1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
     412           2 :                               1.0_dp, subm_tmp1, 'N')
     413             : 
     414             :          !!!!!!! use intermediate matrices to get the error vector !!!!!!!
     415             :          !!!!!!! make sure s_blk_sqrt and its inverse exist (i.e. we use diag algorithm)
     416             :          !CPPrecondition(almo_scf_env%almo_update_algorithm.eq.almo_scf_diag,cp_failure_level,routineP,failure)
     417             :          !! tmp_err = (1-S.T_blk.SigInv.tr(T_blk)).F.T_blk.SigInv
     418             :          !CALL dbcsr_init(matrix_tmp_err)
     419             :          !CALL dbcsr_create(matrix_tmp_err,&
     420             :          !        template=almo_scf_env%matrix_t(ispin))
     421             :          !CALL dbcsr_copy(matrix_tmp_err,&
     422             :          !        matrix_tmp2)
     423             :          !CALL dbcsr_add(matrix_tmp_err,matrix_tmp3,&
     424             :          !        1.0_dp,-1.0_dp)
     425             :          !! err_blk = tmp_err.tr(T_blk)
     426             :          !CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin),&
     427             :          !        almo_scf_env%matrix_s_blk_sqrt(1))
     428             :          !CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp_err,&
     429             :          !        almo_scf_env%matrix_t(ispin),&
     430             :          !        0.0_dp, almo_scf_env%matrix_err_blk(ispin),&
     431             :          !        retain_sparsity=.TRUE.,&
     432             :          !        filter_eps=eps_multiply)
     433             :          !CALL dbcsr_release(matrix_tmp_err)
     434             :          !! bring to the orthogonal basis
     435             :          !! err_blk = (S_blk^-1/2).err_blk.(S_blk^1/2)
     436             :          !CALL dbcsr_init(matrix_tmp_err)
     437             :          !CALL dbcsr_create(matrix_tmp_err,&
     438             :          !        template=almo_scf_env%matrix_err_blk(ispin))
     439             :          !CALL dbcsr_multiply("N", "N", 1.0_dp,&
     440             :          !        almo_scf_env%matrix_err_blk(ispin),&
     441             :          !        almo_scf_env%matrix_s_blk_sqrt(1),&
     442             :          !        0.0_dp, matrix_tmp_err,&
     443             :          !        filter_eps=eps_multiply)
     444             :          !CALL dbcsr_multiply("N", "N", 1.0_dp,&
     445             :          !        almo_scf_env%matrix_s_blk_sqrt_inv(1),&
     446             :          !        matrix_tmp_err,&
     447             :          !        0.0_dp, almo_scf_env%matrix_err_blk(ispin),&
     448             :          !        filter_eps=eps_multiply)
     449             :          !! subtract transpose
     450             :          !CALL dbcsr_transposed(matrix_tmp_err,&
     451             :          !        almo_scf_env%matrix_err_blk(ispin))
     452             :          !CALL dbcsr_add(almo_scf_env%matrix_err_blk(ispin),&
     453             :          !        matrix_tmp_err,&
     454             :          !        1.0_dp,-1.0_dp)
     455             :          !CALL dbcsr_release(matrix_tmp_err)
     456             :          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     457             : 
     458           2 :          CALL release_submatrices(subm_tmp3)
     459           2 :          CALL release_submatrices(subm_tmp2)
     460           2 :          CALL release_submatrices(subm_tmp1)
     461          12 :          DEALLOCATE (subm_tmp3)
     462          12 :          DEALLOCATE (subm_tmp2)
     463          12 :          DEALLOCATE (subm_tmp1)
     464           2 :          CALL dbcsr_release(matrix_tmp3)
     465           2 :          CALL dbcsr_release(matrix_tmp2)
     466           4 :          CALL dbcsr_release(matrix_tmp1)
     467             : 
     468             :       END DO ! spins
     469             : 
     470           2 :       CALL timestop(handle)
     471             : 
     472           4 :    END SUBROUTINE almo_scf_ks_to_ks_xx
     473             : 
     474             : ! **************************************************************************************************
     475             : !> \brief computes the projected KS from the total KS matrix
     476             : !>        also computes the DIIS error vector as a by-product
     477             : !> \param almo_scf_env ...
     478             : !> \par History
     479             : !>       2011.06 created [Rustam Z Khaliullin]
     480             : !> \author Rustam Z Khaliullin
     481             : ! **************************************************************************************************
     482         424 :    SUBROUTINE almo_scf_ks_to_ks_blk(almo_scf_env)
     483             : 
     484             :       TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
     485             : 
     486             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'almo_scf_ks_to_ks_blk'
     487             : 
     488             :       INTEGER                                            :: handle, ispin
     489             :       REAL(KIND=dp)                                      :: eps_multiply
     490             :       TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2, matrix_tmp3, matrix_tmp4, matrix_tmp5, &
     491             :          matrix_tmp6, matrix_tmp7, matrix_tmp8, matrix_tmp9, matrix_tmp_err
     492             : 
     493         424 :       CALL timeset(routineN, handle)
     494             : 
     495         424 :       eps_multiply = almo_scf_env%eps_filter
     496             : 
     497         848 :       DO ispin = 1, almo_scf_env%nspins
     498             : 
     499             :          ! 1. TMP1=KS.T_blk
     500             :          !    Cost: NOn
     501             :          !matrix_tmp1 = create NxO, full
     502             :          CALL dbcsr_create(matrix_tmp1, &
     503         424 :                            template=almo_scf_env%matrix_t(ispin))
     504             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_ks(ispin), &
     505             :                              almo_scf_env%matrix_t_blk(ispin), &
     506             :                              0.0_dp, matrix_tmp1, &
     507         424 :                              filter_eps=eps_multiply)
     508             :          ! 2. TMP2=TMP1.SigInv=KS.T_blk.SigInv
     509             :          !    Cost: NOO
     510             :          !matrix_tmp2 = create NxO, full
     511             :          CALL dbcsr_create(matrix_tmp2, &
     512         424 :                            template=almo_scf_env%matrix_t(ispin))
     513             :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, &
     514             :                              almo_scf_env%matrix_sigma_inv(ispin), &
     515             :                              0.0_dp, matrix_tmp2, &
     516         424 :                              filter_eps=eps_multiply)
     517             : 
     518             :          !!!!!! use intermediate matrices to get the error vector !!!!!!!
     519             :          !CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin),&
     520             :          !        almo_scf_env%matrix_t_blk(ispin))
     521             :          !CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin),&
     522             :          !        matrix_tmp2,&
     523             :          !        keep_sparsity=.TRUE.)
     524             :          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     525             : 
     526             :          ! 3. TMP1=S.T_blk
     527             :          !    Cost: NOn
     528             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s(1), &
     529             :                              almo_scf_env%matrix_t_blk(ispin), &
     530             :                              0.0_dp, matrix_tmp1, &
     531         424 :                              filter_eps=eps_multiply)
     532             : 
     533             :          ! 4. TMP4_blk=TMP2.tr(TMP1)=KS.T_blk.SigInv.tr(T_blk).S
     534             :          !    Cost: NnO
     535             :          !matrix_tmp4 = create NxN, blk
     536             :          CALL dbcsr_create(matrix_tmp4, &
     537             :                            template=almo_scf_env%matrix_s_blk(1), &
     538         424 :                            matrix_type=dbcsr_type_no_symmetry)
     539         424 :          CALL dbcsr_copy(matrix_tmp4, almo_scf_env%matrix_s_blk(1))
     540             :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp2, &
     541             :                              matrix_tmp1, &
     542             :                              0.0_dp, matrix_tmp4, &
     543             :                              retain_sparsity=.TRUE., &
     544         424 :                              filter_eps=eps_multiply)
     545             : 
     546             :          ! 5. KS_blk=KS_blk-TMP4_blk
     547             :          CALL dbcsr_copy(almo_scf_env%matrix_ks_blk(ispin), &
     548         424 :                          almo_scf_env%matrix_ks(ispin), keep_sparsity=.TRUE.)
     549             :          CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), &
     550             :                         matrix_tmp4, &
     551         424 :                         1.0_dp, -1.0_dp)
     552             : 
     553             :          ! 6. TMP5_blk=tr(TMP4_blk)
     554             :          !    KS_blk=KS_blk-tr(TMP4_blk)
     555             :          !matrix_tmp5 = create NxN, blk
     556             :          CALL dbcsr_create(matrix_tmp5, &
     557             :                            template=almo_scf_env%matrix_s_blk(1), &
     558         424 :                            matrix_type=dbcsr_type_no_symmetry)
     559         424 :          CALL dbcsr_transposed(matrix_tmp5, matrix_tmp4)
     560             :          CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp5, &
     561         424 :                         1.0_dp, -1.0_dp)
     562             : 
     563             :          ! 7. TMP3=tr(T_blk).TMP2=tr(T_blk).KS.T_blk.SigInv
     564             :          !    Cost: OOn
     565             :          !matrix_tmp3 = create OxO, full
     566             :          CALL dbcsr_create(matrix_tmp3, &
     567             :                            template=almo_scf_env%matrix_sigma_inv(ispin), &
     568         424 :                            matrix_type=dbcsr_type_no_symmetry)
     569             :          CALL dbcsr_multiply("T", "N", 1.0_dp, &
     570             :                              almo_scf_env%matrix_t_blk(ispin), &
     571             :                              matrix_tmp2, &
     572             :                              0.0_dp, matrix_tmp3, &
     573         424 :                              filter_eps=eps_multiply)
     574             : 
     575             :          ! 8. TMP6=SigInv.TMP3=SigInv.tr(T_blk).KS.T_blk.SigInv
     576             :          !    Cost: OOO
     577             :          !matrix_tmp6 = create OxO, full
     578             :          CALL dbcsr_create(matrix_tmp6, &
     579             :                            template=almo_scf_env%matrix_sigma_inv(ispin), &
     580         424 :                            matrix_type=dbcsr_type_no_symmetry)
     581             :          CALL dbcsr_multiply("N", "N", 1.0_dp, &
     582             :                              almo_scf_env%matrix_sigma_inv(ispin), &
     583             :                              matrix_tmp3, &
     584             :                              0.0_dp, matrix_tmp6, &
     585         424 :                              filter_eps=eps_multiply)
     586             : 
     587             :          ! 9. TMP3=TMP1.TMP6=S.T_blk.SigInv.tr(T_blk).KS.T_blk.SigInv
     588             :          !    Cost: NOO
     589             :          !matrix_tmp3 = re-create NxO, full
     590         424 :          CALL dbcsr_release(matrix_tmp3)
     591             :          CALL dbcsr_create(matrix_tmp3, &
     592         424 :                            template=almo_scf_env%matrix_t(ispin))
     593             :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, &
     594             :                              matrix_tmp6, &
     595             :                              0.0_dp, matrix_tmp3, &
     596         424 :                              filter_eps=eps_multiply)
     597             : 
     598             :          !!!!!! use intermediate matrices to get the error vector !!!!!!!
     599             :          !CALL dbcsr_init(matrix_tmp_err)
     600             :          !CALL dbcsr_create(matrix_tmp_err,&
     601             :          !        template=almo_scf_env%matrix_t_blk(ispin))
     602             :          !CALL dbcsr_copy(matrix_tmp_err,&
     603             :          !        almo_scf_env%matrix_t_blk(ispin))
     604             :          !CALL dbcsr_copy(matrix_tmp_err,matrix_tmp3,&
     605             :          !        keep_sparsity=.TRUE.)
     606             :          !CALL dbcsr_add(almo_scf_env%matrix_err_blk(ispin),matrix_tmp_err,&
     607             :          !        1.0_dp,-1.0_dp)
     608             :          !CALL dbcsr_release(matrix_tmp_err)
     609             :          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     610             : 
     611             :          !!!!!! use intermediate matrices to get the error vector !!!!!!!
     612             :          !!!!!! make sure s_blk_sqrt and its inverse exist (i.e. we use diag algorithm)
     613         424 :          CPASSERT(almo_scf_env%almo_update_algorithm .EQ. almo_scf_diag)
     614             :          ! tmp_err = (1-S.T_blk.SigInv.tr(T_blk)).F.T_blk.SigInv
     615             :          CALL dbcsr_create(matrix_tmp_err, &
     616         424 :                            template=almo_scf_env%matrix_t_blk(ispin))
     617             :          CALL dbcsr_copy(matrix_tmp_err, &
     618         424 :                          matrix_tmp2)
     619             :          CALL dbcsr_add(matrix_tmp_err, matrix_tmp3, &
     620         424 :                         1.0_dp, -1.0_dp)
     621             :          ! err_blk = tmp_err.tr(T_blk)
     622             :          CALL dbcsr_copy(almo_scf_env%matrix_err_blk(ispin), &
     623         424 :                          almo_scf_env%matrix_s_blk_sqrt(1))
     624             :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp_err, &
     625             :                              almo_scf_env%matrix_t_blk(ispin), &
     626             :                              0.0_dp, almo_scf_env%matrix_err_blk(ispin), &
     627             :                              retain_sparsity=.TRUE., &
     628         424 :                              filter_eps=eps_multiply)
     629         424 :          CALL dbcsr_release(matrix_tmp_err)
     630             :          ! bring to the orthogonal basis
     631             :          ! err_blk = (S_blk^-1/2).err_blk.(S_blk^1/2)
     632             :          CALL dbcsr_create(matrix_tmp_err, &
     633         424 :                            template=almo_scf_env%matrix_err_blk(ispin))
     634             :          CALL dbcsr_multiply("N", "N", 1.0_dp, &
     635             :                              almo_scf_env%matrix_err_blk(ispin), &
     636             :                              almo_scf_env%matrix_s_blk_sqrt(1), &
     637             :                              0.0_dp, matrix_tmp_err, &
     638         424 :                              filter_eps=eps_multiply)
     639             :          CALL dbcsr_multiply("N", "N", 1.0_dp, &
     640             :                              almo_scf_env%matrix_s_blk_sqrt_inv(1), &
     641             :                              matrix_tmp_err, &
     642             :                              0.0_dp, almo_scf_env%matrix_err_blk(ispin), &
     643         424 :                              filter_eps=eps_multiply)
     644             : 
     645             :          ! subtract transpose
     646             :          CALL dbcsr_transposed(matrix_tmp_err, &
     647         424 :                                almo_scf_env%matrix_err_blk(ispin))
     648             :          CALL dbcsr_add(almo_scf_env%matrix_err_blk(ispin), &
     649             :                         matrix_tmp_err, &
     650         424 :                         1.0_dp, -1.0_dp)
     651         424 :          CALL dbcsr_release(matrix_tmp_err)
     652             :          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     653             : 
     654             :          ! later we will need only the blk version of TMP6
     655             :          ! create it here and release TMP6
     656             :          !matrix_tmp9 = create OxO, blk
     657             :          !matrix_tmp9 = copy data from matrix_tmp6, retain sparsity
     658             :          !matrix_tmp6 = release
     659             :          CALL dbcsr_create(matrix_tmp9, &
     660             :                            template=almo_scf_env%matrix_sigma_blk(ispin), &
     661         424 :                            matrix_type=dbcsr_type_no_symmetry)
     662         424 :          CALL dbcsr_copy(matrix_tmp9, almo_scf_env%matrix_sigma_blk(ispin))
     663         424 :          CALL dbcsr_copy(matrix_tmp9, matrix_tmp6, keep_sparsity=.TRUE.)
     664         424 :          CALL dbcsr_release(matrix_tmp6)
     665             : 
     666             :          !10. KS_blk=KS_blk+TMP3.tr(TMP1)=
     667             :          !          =KS_blk+S.T_blk.SigInv.tr.(T_blk).KS.T_blk.SigInv.tr(T_blk).S
     668             :          !    Cost: NnO
     669             :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp3, &
     670             :                              matrix_tmp1, &
     671             :                              1.0_dp, almo_scf_env%matrix_ks_blk(ispin), &
     672             :                              retain_sparsity=.TRUE., &
     673         424 :                              filter_eps=eps_multiply)
     674             : 
     675             :          ! 11. TMP4_blk=TMP7_blk.tr(TMP8_blk)
     676             :          !    Cost: Nnn
     677             :          !matrix_tmp7 = create NxO, blk
     678             :          !matrix_tmp7 = copy data from matrix_tmp3, retain sparsity
     679             :          !matrix_tmp3 = release
     680             :          !matrix_tmp8 = create NxO, blk
     681             :          !matrix_tmp8 = copy data from matrix_tmp1, retain sparsity
     682             :          !matrix_tmp1 = release
     683             :          CALL dbcsr_create(matrix_tmp7, &
     684         424 :                            template=almo_scf_env%matrix_t_blk(ispin))
     685             :          ! transfer only the ALMO blocks from tmp3 into tmp7:
     686             :          ! first, copy t_blk into tmp7 to transfer the blk structure,
     687             :          ! then copy tmp3 into tmp7 with retain_sparsity
     688         424 :          CALL dbcsr_copy(matrix_tmp7, almo_scf_env%matrix_t_blk(ispin))
     689         424 :          CALL dbcsr_copy(matrix_tmp7, matrix_tmp3, keep_sparsity=.TRUE.)
     690         424 :          CALL dbcsr_release(matrix_tmp3)
     691             :          ! do the same for tmp1->tmp8
     692             :          CALL dbcsr_create(matrix_tmp8, &
     693         424 :                            template=almo_scf_env%matrix_t_blk(ispin))
     694         424 :          CALL dbcsr_copy(matrix_tmp8, almo_scf_env%matrix_t_blk(ispin))
     695         424 :          CALL dbcsr_copy(matrix_tmp8, matrix_tmp1, keep_sparsity=.TRUE.)
     696         424 :          CALL dbcsr_release(matrix_tmp1)
     697             :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp7, &
     698             :                              matrix_tmp8, &
     699             :                              0.0_dp, matrix_tmp4, &
     700             :                              filter_eps=eps_multiply, &
     701         424 :                              retain_sparsity=.TRUE.)
     702             : 
     703             :          ! 12. KS_blk=KS_blk-TMP4_blk
     704             :          CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp4, &
     705         424 :                         1.0_dp, -1.0_dp)
     706             : 
     707             :          ! 13. TMP5_blk=tr(TMP5_blk)
     708             :          !     KS_blk=KS_blk-tr(TMP4_blk)
     709         424 :          CALL dbcsr_transposed(matrix_tmp5, matrix_tmp4)
     710             :          CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp5, &
     711         424 :                         1.0_dp, -1.0_dp)
     712             : 
     713             :          ! 14. TMP4_blk=TMP7_blk.tr(TMP8_blk)
     714             :          !     Cost: Nnn
     715         424 :          CALL dbcsr_copy(matrix_tmp7, matrix_tmp2, keep_sparsity=.TRUE.)
     716         424 :          CALL dbcsr_release(matrix_tmp2)
     717             :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp7, &
     718             :                              matrix_tmp8, &
     719             :                              0.0_dp, matrix_tmp4, &
     720             :                              retain_sparsity=.TRUE., &
     721         424 :                              filter_eps=eps_multiply)
     722             :          ! 15. KS_blk=KS_blk+TMP4_blk
     723             :          CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp4, &
     724         424 :                         1.0_dp, 1.0_dp)
     725             : 
     726             :          ! 16. KS_blk=KS_blk+tr(TMP4_blk)
     727         424 :          CALL dbcsr_transposed(matrix_tmp5, matrix_tmp4)
     728         424 :          CALL dbcsr_release(matrix_tmp4)
     729             :          CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), matrix_tmp5, &
     730         424 :                         1.0_dp, 1.0_dp)
     731         424 :          CALL dbcsr_release(matrix_tmp5)
     732             : 
     733             :          ! 17. TMP10_blk=TMP8_blk.TMP9_blk
     734             :          !    Cost: Noo
     735             :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp8, &
     736             :                              matrix_tmp9, &
     737             :                              0.0_dp, matrix_tmp7, &
     738             :                              retain_sparsity=.TRUE., &
     739         424 :                              filter_eps=eps_multiply)
     740         424 :          CALL dbcsr_release(matrix_tmp9)
     741             : 
     742             :          ! 18. KS_blk=TMP7_blk.tr(TMP8_blk)
     743             :          !    Cost: Nno
     744             :          CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp7, &
     745             :                              matrix_tmp8, &
     746             :                              1.0_dp, almo_scf_env%matrix_ks_blk(ispin), &
     747             :                              retain_sparsity=.TRUE., &
     748         424 :                              filter_eps=eps_multiply)
     749         424 :          CALL dbcsr_release(matrix_tmp7)
     750         848 :          CALL dbcsr_release(matrix_tmp8)
     751             : 
     752             :       END DO ! spins
     753             : 
     754         424 :       CALL timestop(handle)
     755             : 
     756         424 :    END SUBROUTINE almo_scf_ks_to_ks_blk
     757             : 
     758             : ! **************************************************************************************************
     759             : !> \brief ALMOs by diagonalizing the KS domain submatrices
     760             : !>        computes both the occupied and virtual orbitals
     761             : !> \param almo_scf_env ...
     762             : !> \par History
     763             : !>       2013.03 created [Rustam Z Khaliullin]
     764             : !>       2018.09 smearing support [Ruben Staub]
     765             : !> \author Rustam Z Khaliullin
     766             : ! **************************************************************************************************
     767           2 :    SUBROUTINE almo_scf_ks_xx_to_tv_xx(almo_scf_env)
     768             : 
     769             :       TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
     770             : 
     771             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'almo_scf_ks_xx_to_tv_xx'
     772             : 
     773             :       INTEGER                                            :: handle, iblock_size, idomain, info, &
     774             :                                                             ispin, lwork, ndomains
     775           2 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, work
     776           2 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: data_copy
     777             :       TYPE(domain_submatrix_type), ALLOCATABLE, &
     778           2 :          DIMENSION(:)                                    :: subm_ks_xx_orthog, subm_t, subm_tmp
     779             : 
     780           2 :       CALL timeset(routineN, handle)
     781             : 
     782           2 :       IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular .AND. &
     783             :           almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
     784           0 :          CPABORT("a domain must be located entirely on a CPU")
     785             :       END IF
     786             : 
     787           2 :       ndomains = almo_scf_env%ndomains
     788          16 :       ALLOCATE (subm_tmp(ndomains))
     789          16 :       ALLOCATE (subm_ks_xx_orthog(ndomains))
     790          16 :       ALLOCATE (subm_t(ndomains))
     791             : 
     792           4 :       DO ispin = 1, almo_scf_env%nspins
     793             : 
     794           2 :          CALL init_submatrices(subm_tmp)
     795           2 :          CALL init_submatrices(subm_ks_xx_orthog)
     796             : 
     797             :          ! TRY: project out T0-occupied space for each domain
     798             :          ! F=(1-R_du).F.(1-tr(R_du))
     799             :          !CALL copy_submatrices(almo_scf_env%domain_ks_xx(:,ispin),&
     800             :          !        subm_ks_xx_orthog,copy_data=.TRUE.)
     801             :          !CALL multiply_submatrices('N','N',1.0_dp,&
     802             :          !        almo_scf_env%domain_r_down_up(:,ispin),&
     803             :          !        almo_scf_env%domain_ks_xx(:,ispin),0.0_dp,subm_tmp)
     804             :          !CALL add_submatrices(1.0_dp,subm_ks_xx_orthog,-1.0_dp,subm_tmp,'N')
     805             :          !CALL add_submatrices(1.0_dp,subm_ks_xx_orthog,-1.0_dp,subm_tmp,'T')
     806             :          !!CALL multiply_submatrices('N','T',1.0_dp,subm_tmp,&
     807             :          !!        almo_scf_env%domain_r_down_up(:,ispin),&
     808             :          !!        1.0_dp,subm_ks_xx_orthog)
     809             : 
     810             :          ! convert blocks to the orthogonal basis set
     811             :          ! TRY: replace one multiply
     812             :          !CALL multiply_submatrices('N','N',1.0_dp,subm_ks_xx_orthog,&
     813             :          !        almo_scf_env%domain_s_sqrt_inv(:,ispin),0.0_dp,subm_tmp)
     814             :          CALL multiply_submatrices('N', 'N', 1.0_dp, almo_scf_env%domain_ks_xx(:, ispin), &
     815           2 :                                    almo_scf_env%domain_s_sqrt_inv(:, ispin), 0.0_dp, subm_tmp)
     816             :          CALL multiply_submatrices('N', 'N', 1.0_dp, almo_scf_env%domain_s_sqrt_inv(:, ispin), &
     817           2 :                                    subm_tmp, 0.0_dp, subm_ks_xx_orthog)
     818           2 :          CALL release_submatrices(subm_tmp)
     819             : 
     820             :          ! create temporary matrices for occupied and virtual orbitals
     821             :          ! represented in the orthogonalized basis set
     822           2 :          CALL init_submatrices(subm_t)
     823             : 
     824             :          ! loop over domains - perform diagonalization
     825          12 :          DO idomain = 1, ndomains
     826             : 
     827             :             ! check if the submatrix exists
     828          12 :             IF (subm_ks_xx_orthog(idomain)%domain .GT. 0) THEN
     829             : 
     830           5 :                iblock_size = subm_ks_xx_orthog(idomain)%nrows
     831             : 
     832             :                ! Prepare data
     833          15 :                ALLOCATE (eigenvalues(iblock_size))
     834          20 :                ALLOCATE (data_copy(iblock_size, iblock_size))
     835       31607 :                data_copy(:, :) = subm_ks_xx_orthog(idomain)%mdata(:, :)
     836             : 
     837             :                ! Query the optimal workspace for dsyev
     838           5 :                LWORK = -1
     839           5 :                ALLOCATE (WORK(MAX(1, LWORK)))
     840           5 :                CALL DSYEV('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
     841           5 :                LWORK = INT(WORK(1))
     842           5 :                DEALLOCATE (WORK)
     843             : 
     844             :                ! Allocate the workspace and solve the eigenproblem
     845          15 :                ALLOCATE (WORK(MAX(1, LWORK)))
     846           5 :                CALL DSYEV('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
     847           5 :                IF (INFO .NE. 0) THEN
     848           0 :                   CPABORT("DSYEV failed")
     849             :                END IF
     850             : 
     851             :                ! Copy occupied eigenvectors
     852           5 :                IF (almo_scf_env%domain_t(idomain, ispin)%ncols .NE. &
     853             :                    almo_scf_env%nocc_of_domain(idomain, ispin)) THEN
     854           0 :                   CPABORT("wrong domain structure")
     855             :                END IF
     856             :                CALL copy_submatrices(almo_scf_env%domain_t(idomain, ispin), &
     857           5 :                                      subm_t(idomain), .FALSE.)
     858             :                CALL copy_submatrix_data(data_copy(:, 1:almo_scf_env%nocc_of_domain(idomain, ispin)), &
     859           5 :                                         subm_t(idomain))
     860             :                !! Copy occupied eigenvalues if smearing requested
     861           5 :                IF (almo_scf_env%smear) THEN
     862             :                   almo_scf_env%mo_energies(1 + SUM(almo_scf_env%nocc_of_domain(:idomain - 1, ispin)) &
     863             :                                            :SUM(almo_scf_env%nocc_of_domain(:idomain, ispin)), ispin) &
     864           0 :                      = eigenvalues(1:almo_scf_env%nocc_of_domain(idomain, ispin))
     865             :                END IF
     866             : 
     867           5 :                DEALLOCATE (WORK)
     868           5 :                DEALLOCATE (data_copy)
     869           5 :                DEALLOCATE (eigenvalues)
     870             : 
     871             :             END IF ! submatrix for the domain exists
     872             : 
     873             :          END DO ! loop over domains
     874             : 
     875           2 :          CALL release_submatrices(subm_ks_xx_orthog)
     876             : 
     877             :          ! convert orbitals to the AO basis set (from orthogonalized AOs)
     878             :          CALL multiply_submatrices('N', 'N', 1.0_dp, almo_scf_env%domain_s_sqrt_inv(:, ispin), &
     879           2 :                                    subm_t, 0.0_dp, almo_scf_env%domain_t(:, ispin))
     880           2 :          CALL release_submatrices(subm_t)
     881             : 
     882             :          ! convert domain orbitals to a dbcsr matrix
     883             :          CALL construct_dbcsr_from_submatrices( &
     884             :             almo_scf_env%matrix_t(ispin), &
     885             :             almo_scf_env%domain_t(:, ispin), &
     886           2 :             almo_scf_env%quench_t(ispin))
     887             :          CALL dbcsr_filter(almo_scf_env%matrix_t(ispin), &
     888           4 :                            almo_scf_env%eps_filter)
     889             : 
     890             :          ! TRY: add T0 component
     891             :          !!CALL dbcsr_add(almo_scf_env%matrix_t(ispin),&
     892             :          !!        almo_scf_env%matrix_t_blk(ispin),1.0_dp,1.0_dp)
     893             : 
     894             :       END DO ! spins
     895             : 
     896          12 :       DEALLOCATE (subm_tmp)
     897          12 :       DEALLOCATE (subm_ks_xx_orthog)
     898          12 :       DEALLOCATE (subm_t)
     899             : 
     900           2 :       CALL timestop(handle)
     901             : 
     902           4 :    END SUBROUTINE almo_scf_ks_xx_to_tv_xx
     903             : 
     904             : ! **************************************************************************************************
     905             : !> \brief computes ALMOs by diagonalizing the projected blocked KS matrix
     906             : !>        uses the diagonalization code for blocks
     907             : !>        computes both the occupied and virtual orbitals
     908             : !> \param almo_scf_env ...
     909             : !> \par History
     910             : !>       2011.07 created [Rustam Z Khaliullin]
     911             : !>       2018.09 smearing support [Ruben Staub]
     912             : !> \author Rustam Z Khaliullin
     913             : ! **************************************************************************************************
     914         348 :    SUBROUTINE almo_scf_ks_blk_to_tv_blk(almo_scf_env)
     915             : 
     916             :       TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
     917             : 
     918             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'almo_scf_ks_blk_to_tv_blk'
     919             : 
     920             :       INTEGER                                            :: handle, iblock_col, iblock_row, &
     921             :                                                             iblock_size, info, ispin, lwork, &
     922             :                                                             nocc_of_block, nvirt_of_block, orbital
     923             :       LOGICAL                                            :: block_needed
     924         348 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, work
     925         348 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: data_copy
     926         348 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: data_p, p_new_block
     927             :       TYPE(dbcsr_iterator_type)                          :: iter
     928             :       TYPE(dbcsr_type)                                   :: matrix_ks_blk_orthog, &
     929             :                                                             matrix_t_blk_orthog, matrix_tmp, &
     930             :                                                             matrix_v_blk_orthog
     931             : 
     932         348 :       CALL timeset(routineN, handle)
     933             : 
     934         348 :       IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular .AND. &
     935             :           almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
     936           0 :          CPABORT("a domain must be located entirely on a CPU")
     937             :       END IF
     938             : 
     939         696 :       DO ispin = 1, almo_scf_env%nspins
     940             : 
     941             :          CALL dbcsr_create(matrix_tmp, template=almo_scf_env%matrix_ks_blk(ispin), &
     942         348 :                            matrix_type=dbcsr_type_no_symmetry)
     943             :          CALL dbcsr_create(matrix_ks_blk_orthog, template=almo_scf_env%matrix_ks_blk(ispin), &
     944         348 :                            matrix_type=dbcsr_type_no_symmetry)
     945             : 
     946             :          ! convert blocks to the orthogonal basis set
     947             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_ks_blk(ispin), &
     948             :                              almo_scf_env%matrix_s_blk_sqrt_inv(1), 0.0_dp, matrix_tmp, &
     949         348 :                              filter_eps=almo_scf_env%eps_filter)
     950             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk_sqrt_inv(1), &
     951             :                              matrix_tmp, 0.0_dp, matrix_ks_blk_orthog, &
     952         348 :                              filter_eps=almo_scf_env%eps_filter)
     953             : 
     954         348 :          CALL dbcsr_release(matrix_tmp)
     955             : 
     956             :          ! create temporary matrices for occupied and virtual orbitals
     957             :          ! represented in the orthogonalized AOs basis set
     958         348 :          CALL dbcsr_create(matrix_t_blk_orthog, template=almo_scf_env%matrix_t_blk(ispin))
     959         348 :          CALL dbcsr_create(matrix_v_blk_orthog, template=almo_scf_env%matrix_v_full_blk(ispin))
     960         348 :          CALL dbcsr_work_create(matrix_t_blk_orthog, work_mutable=.TRUE.)
     961         348 :          CALL dbcsr_work_create(matrix_v_blk_orthog, work_mutable=.TRUE.)
     962             : 
     963         348 :          CALL dbcsr_work_create(almo_scf_env%matrix_eoo(ispin), work_mutable=.TRUE.)
     964         348 :          CALL dbcsr_work_create(almo_scf_env%matrix_evv_full(ispin), work_mutable=.TRUE.)
     965             : 
     966         348 :          CALL dbcsr_iterator_start(iter, matrix_ks_blk_orthog)
     967             : 
     968        1624 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     969        1276 :             CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
     970             : 
     971        1276 :             IF (iblock_row .NE. iblock_col) THEN
     972           0 :                CPABORT("off-diagonal block found")
     973             :             END IF
     974             : 
     975        1276 :             block_needed = .TRUE.
     976        1276 :             IF (almo_scf_env%nocc_of_domain(iblock_col, ispin) .EQ. 0 .AND. &
     977             :                 almo_scf_env%nvirt_of_domain(iblock_col, ispin) .EQ. 0) THEN
     978             :                block_needed = .FALSE.
     979             :             END IF
     980             : 
     981         348 :             IF (block_needed) THEN
     982             : 
     983             :                ! Prepare data
     984        3828 :                ALLOCATE (eigenvalues(iblock_size))
     985        5104 :                ALLOCATE (data_copy(iblock_size, iblock_size))
     986      382688 :                data_copy(:, :) = data_p(:, :)
     987             : 
     988             :                ! Query the optimal workspace for dsyev
     989        1276 :                LWORK = -1
     990        1276 :                ALLOCATE (WORK(MAX(1, LWORK)))
     991        1276 :                CALL DSYEV('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
     992        1276 :                LWORK = INT(WORK(1))
     993        1276 :                DEALLOCATE (WORK)
     994             : 
     995             :                ! Allocate the workspace and solve the eigenproblem
     996        3828 :                ALLOCATE (WORK(MAX(1, LWORK)))
     997        1276 :                CALL DSYEV('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
     998        1276 :                IF (INFO .NE. 0) THEN
     999           0 :                   CPABORT("DSYEV failed")
    1000             :                END IF
    1001             : 
    1002             :                !!! RZK-warning                                               !!!
    1003             :                !!! IT IS EXTREMELY IMPORTANT THAT THE DIAGONAL BLOCKS OF THE !!!
    1004             :                !!! FOLLOWING MATRICES ARE LOCATED ON THE SAME NODES WITH     !!!
    1005             :                !!! THE CORRESPONDING DIAGONAL BLOCKS OF THE FOCK MATRIX:     !!!
    1006             :                !!! T, V, E_o, E_v
    1007             : 
    1008             :                ! copy eigenvectors into two dbcsr matrices - occupied and virtuals
    1009        1276 :                nocc_of_block = almo_scf_env%nocc_of_domain(iblock_col, ispin)
    1010        1276 :                IF (nocc_of_block .GT. 0) THEN
    1011        1252 :                   NULLIFY (p_new_block)
    1012        1252 :                   CALL dbcsr_reserve_block2d(matrix_t_blk_orthog, iblock_row, iblock_col, p_new_block)
    1013        1252 :                   CPASSERT(ASSOCIATED(p_new_block))
    1014       90468 :                   p_new_block(:, :) = data_copy(:, 1:nocc_of_block)
    1015             :                   ! copy eigenvalues into diagonal dbcsr matrix - Eoo
    1016        1252 :                   NULLIFY (p_new_block)
    1017        1252 :                   CALL dbcsr_reserve_block2d(almo_scf_env%matrix_eoo(ispin), iblock_row, iblock_col, p_new_block)
    1018        1252 :                   CPASSERT(ASSOCIATED(p_new_block))
    1019       32156 :                   p_new_block(:, :) = 0.0_dp
    1020        6358 :                   DO orbital = 1, nocc_of_block
    1021        6358 :                      p_new_block(orbital, orbital) = eigenvalues(orbital)
    1022             :                   END DO
    1023             :                   !! Retrieve occupied MOs energies for smearing purpose, if requested
    1024             :                   !! RS-WARNING: Hack to retrieve the occupied energies, since matrix_eoo seems to be ill-defined
    1025             :                   !!             for multiprocessing (any idea for fix?)
    1026             :                   !! RS-WARNING: This section is not suitable for parallel run !!!
    1027             :                   !!             (but usually fails less than retrieving the diagonal of matrix_eoo)
    1028             :                   !! RS-WARNING: This method will likely keep the energies of the initial guess if run in parallel
    1029             :                   !!             (which is still a reasonable smearing in most cases...)
    1030        1252 :                   IF (almo_scf_env%smear) THEN
    1031         292 :                      DO orbital = 1, nocc_of_block
    1032             :                         almo_scf_env%mo_energies(SUM(almo_scf_env%nocc_of_domain(:iblock_row - 1, ispin)) + orbital, &
    1033         348 :                                                  ispin) = eigenvalues(orbital)
    1034             :                      END DO
    1035             :                   END IF
    1036             :                END IF
    1037             : 
    1038             :                ! now virtuals
    1039        1276 :                nvirt_of_block = almo_scf_env%nvirt_of_domain(iblock_col, ispin)
    1040        1276 :                IF (nvirt_of_block .GT. 0) THEN
    1041        1276 :                   NULLIFY (p_new_block)
    1042        1276 :                   CALL dbcsr_reserve_block2d(matrix_v_blk_orthog, iblock_row, iblock_col, p_new_block)
    1043        1276 :                   CPASSERT(ASSOCIATED(p_new_block))
    1044      293472 :                   p_new_block(:, :) = data_copy(:, (nocc_of_block + 1):(nocc_of_block + nvirt_of_block))
    1045             :                   ! virtual energies
    1046        1276 :                   NULLIFY (p_new_block)
    1047        1276 :                   CALL dbcsr_reserve_block2d(almo_scf_env%matrix_evv_full(ispin), iblock_row, iblock_col, p_new_block)
    1048        1276 :                   CPASSERT(ASSOCIATED(p_new_block))
    1049      235160 :                   p_new_block(:, :) = 0.0_dp
    1050       13896 :                   DO orbital = 1, nvirt_of_block
    1051       13896 :                      p_new_block(orbital, orbital) = eigenvalues(nocc_of_block + orbital)
    1052             :                   END DO
    1053             :                END IF
    1054             : 
    1055        1276 :                DEALLOCATE (WORK)
    1056        1276 :                DEALLOCATE (data_copy)
    1057        1276 :                DEALLOCATE (eigenvalues)
    1058             : 
    1059             :             END IF
    1060             : 
    1061             :          END DO
    1062         348 :          CALL dbcsr_iterator_stop(iter)
    1063             : 
    1064         348 :          CALL dbcsr_finalize(matrix_t_blk_orthog)
    1065         348 :          CALL dbcsr_finalize(matrix_v_blk_orthog)
    1066         348 :          CALL dbcsr_finalize(almo_scf_env%matrix_eoo(ispin))
    1067         348 :          CALL dbcsr_finalize(almo_scf_env%matrix_evv_full(ispin))
    1068             : 
    1069             :          !! RS-WARNING: When matrix_eoo will be well-defined with multiprocessing,
    1070             :          !!             the following should be the preferred way to retrieve the occupied energies:
    1071             :          !! Retrieve occupied MOs energies for smearing purpose, if requested
    1072             :          !! IF (almo_scf_env%smear) THEN
    1073             :          !!    CALL dbcsr_get_diag(almo_scf_env%matrix_eoo(ispin), almo_scf_env%mo_energies(:, ispin))
    1074             :          !! END IF
    1075             : 
    1076         348 :          CALL dbcsr_filter(matrix_t_blk_orthog, almo_scf_env%eps_filter)
    1077         348 :          CALL dbcsr_filter(matrix_v_blk_orthog, almo_scf_env%eps_filter)
    1078             : 
    1079         348 :          CALL dbcsr_release(matrix_ks_blk_orthog)
    1080             : 
    1081             :          ! convert orbitals to the AO basis set (from orthogonalized AOs)
    1082             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk_sqrt_inv(1), &
    1083             :                              matrix_t_blk_orthog, 0.0_dp, almo_scf_env%matrix_t_blk(ispin), &
    1084         348 :                              filter_eps=almo_scf_env%eps_filter)
    1085             :          CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk_sqrt_inv(1), &
    1086             :                              matrix_v_blk_orthog, 0.0_dp, almo_scf_env%matrix_v_full_blk(ispin), &
    1087         348 :                              filter_eps=almo_scf_env%eps_filter)
    1088             : 
    1089         348 :          CALL dbcsr_release(matrix_t_blk_orthog)
    1090         696 :          CALL dbcsr_release(matrix_v_blk_orthog)
    1091             : 
    1092             :       END DO ! spins
    1093             : 
    1094         348 :       CALL timestop(handle)
    1095             : 
    1096         696 :    END SUBROUTINE almo_scf_ks_blk_to_tv_blk
    1097             : 
    1098             : ! **************************************************************************************************
    1099             : !> \brief inverts block-diagonal blocks of a dbcsr_matrix
    1100             : !> \param matrix_in ...
    1101             : !> \param matrix_out ...
    1102             : !> \param nocc ...
    1103             : !> \par History
    1104             : !>       2012.05 created [Rustam Z Khaliullin]
    1105             : !> \author Rustam Z Khaliullin
    1106             : ! **************************************************************************************************
    1107          82 :    SUBROUTINE pseudo_invert_diagonal_blk(matrix_in, matrix_out, nocc)
    1108             : 
    1109             :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_in
    1110             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_out
    1111             :       INTEGER, DIMENSION(:)                              :: nocc
    1112             : 
    1113             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'pseudo_invert_diagonal_blk'
    1114             : 
    1115             :       INTEGER                                            :: handle, iblock_col, iblock_row, &
    1116             :                                                             iblock_size, methodID
    1117          82 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: data_copy
    1118          82 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: data_p, p_new_block
    1119             :       TYPE(dbcsr_iterator_type)                          :: iter
    1120             : 
    1121          82 :       CALL timeset(routineN, handle)
    1122             : 
    1123          82 :       CALL dbcsr_create(matrix_out, template=matrix_in)
    1124          82 :       CALL dbcsr_work_create(matrix_out, work_mutable=.TRUE.)
    1125             : 
    1126          82 :       CALL dbcsr_iterator_start(iter, matrix_in)
    1127             : 
    1128         423 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
    1129             : 
    1130         341 :          CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
    1131             : 
    1132         423 :          IF (iblock_row == iblock_col) THEN
    1133             : 
    1134             :             ! Prepare data
    1135        1276 :             ALLOCATE (data_copy(iblock_size, iblock_size))
    1136             :             !data_copy(:,:)=data_p(:,:)
    1137             : 
    1138             :             ! 0. Cholesky factorization
    1139             :             ! 1. Diagonalization
    1140         319 :             methodID = 1
    1141             :             CALL pseudo_invert_matrix(data_p, data_copy, iblock_size, &
    1142             :                                       methodID, &
    1143             :                                       range1=nocc(iblock_row), range2=nocc(iblock_row), &
    1144             :                                       !range1_thr,range2_thr,&
    1145         319 :                                       shift=1.0E-5_dp)
    1146             :             !!! IT IS EXTREMELY IMPORTANT THAT THE BLOCKS OF THE "OUT"  !!!
    1147             :             !!! MATRIX ARE DISTRIBUTED AS THE BLOCKS OF THE "IN" MATRIX !!!
    1148             : 
    1149         319 :             NULLIFY (p_new_block)
    1150         319 :             CALL dbcsr_reserve_block2d(matrix_out, iblock_row, iblock_col, p_new_block)
    1151         319 :             CPASSERT(ASSOCIATED(p_new_block))
    1152      221823 :             p_new_block(:, :) = data_copy(:, :)
    1153             : 
    1154         319 :             DEALLOCATE (data_copy)
    1155             : 
    1156             :          END IF
    1157             : 
    1158             :       END DO
    1159          82 :       CALL dbcsr_iterator_stop(iter)
    1160             : 
    1161          82 :       CALL dbcsr_finalize(matrix_out)
    1162             : 
    1163          82 :       CALL timestop(handle)
    1164             : 
    1165         164 :    END SUBROUTINE pseudo_invert_diagonal_blk
    1166             : 
    1167             : ! **************************************************************************************************
    1168             : !> \brief computes occupied ALMOs from the superimposed atomic density blocks
    1169             : !> \param almo_scf_env ...
    1170             : !> \param ionic ...
    1171             : !> \par History
    1172             : !>       2011.06 created [Rustam Z Khaliullin]
    1173             : !> \author Rustam Z Khaliullin
    1174             : ! **************************************************************************************************
    1175          60 :    SUBROUTINE almo_scf_p_blk_to_t_blk(almo_scf_env, ionic)
    1176             : 
    1177             :       TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
    1178             :       LOGICAL, INTENT(IN)                                :: ionic
    1179             : 
    1180             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'almo_scf_p_blk_to_t_blk'
    1181             : 
    1182             :       INTEGER                                            :: handle, iblock_col, iblock_row, &
    1183             :                                                             iblock_size, info, ispin, lwork, &
    1184             :                                                             nocc_of_block, unit_nr
    1185             :       LOGICAL                                            :: block_needed
    1186          60 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, work
    1187          60 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: data_copy
    1188          60 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: data_p, p_new_block
    1189             :       TYPE(cp_logger_type), POINTER                      :: logger
    1190             :       TYPE(dbcsr_iterator_type)                          :: iter
    1191             :       TYPE(dbcsr_type)                                   :: matrix_t_blk_tmp
    1192             : 
    1193          60 :       CALL timeset(routineN, handle)
    1194             : 
    1195             :       ! get a useful unit_nr
    1196          60 :       logger => cp_get_default_logger()
    1197          60 :       IF (logger%para_env%is_source()) THEN
    1198          30 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    1199             :       ELSE
    1200             :          unit_nr = -1
    1201             :       END IF
    1202             : 
    1203         120 :       DO ispin = 1, almo_scf_env%nspins
    1204             : 
    1205         120 :          IF (ionic) THEN
    1206             : 
    1207             :             ! create a temporary matrix to keep the eigenvectors
    1208             :             CALL dbcsr_create(matrix_t_blk_tmp, &
    1209           0 :                               template=almo_scf_env%matrix_t_blk(ispin))
    1210             :             CALL dbcsr_work_create(matrix_t_blk_tmp, &
    1211           0 :                                    work_mutable=.TRUE.)
    1212             : 
    1213           0 :             CALL dbcsr_iterator_start(iter, almo_scf_env%matrix_p_blk(ispin))
    1214           0 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
    1215           0 :                CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
    1216             : 
    1217           0 :                block_needed = .FALSE.
    1218             : 
    1219           0 :                IF (iblock_row == iblock_col) THEN
    1220             :                   block_needed = .TRUE.
    1221             :                END IF
    1222             : 
    1223             :                IF (.NOT. block_needed) THEN
    1224           0 :                   CPABORT("off-diag block found")
    1225             :                END IF
    1226             : 
    1227             :                ! Prepare data
    1228           0 :                ALLOCATE (eigenvalues(iblock_size))
    1229           0 :                ALLOCATE (data_copy(iblock_size, iblock_size))
    1230           0 :                data_copy(:, :) = data_p(:, :)
    1231             : 
    1232             :                ! Query the optimal workspace for dsyev
    1233           0 :                LWORK = -1
    1234           0 :                ALLOCATE (WORK(MAX(1, LWORK)))
    1235             :                CALL DSYEV('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, &
    1236           0 :                           WORK, LWORK, INFO)
    1237           0 :                LWORK = INT(WORK(1))
    1238           0 :                DEALLOCATE (WORK)
    1239             : 
    1240             :                ! Allocate the workspace and solve the eigenproblem
    1241           0 :                ALLOCATE (WORK(MAX(1, LWORK)))
    1242           0 :                CALL DSYEV('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
    1243           0 :                IF (INFO .NE. 0) THEN
    1244           0 :                   IF (unit_nr > 0) THEN
    1245           0 :                      WRITE (unit_nr, *) 'BLOCK  = ', iblock_row
    1246           0 :                      WRITE (unit_nr, *) 'INFO   =', INFO
    1247           0 :                      WRITE (unit_nr, *) data_p(:, :)
    1248             :                   END IF
    1249           0 :                   CPABORT("DSYEV failed")
    1250             :                END IF
    1251             : 
    1252             :                !!! IT IS EXTREMELY IMPORTANT THAT THE DIAGONAL BLOCKS OF THE !!!
    1253             :                !!! P AND T MATRICES ARE LOCATED ON THE SAME NODES            !!!
    1254             : 
    1255             :                ! copy eigenvectors into two dbcsr matrices - occupied and virtuals
    1256           0 :                NULLIFY (p_new_block)
    1257             :                CALL dbcsr_reserve_block2d(matrix_t_blk_tmp, &
    1258           0 :                                           iblock_row, iblock_col, p_new_block)
    1259           0 :                nocc_of_block = SIZE(p_new_block, 2)
    1260           0 :                CPASSERT(ASSOCIATED(p_new_block))
    1261           0 :                CPASSERT(nocc_of_block .GT. 0)
    1262           0 :                p_new_block(:, :) = data_copy(:, iblock_size - nocc_of_block + 1:)
    1263             : 
    1264           0 :                DEALLOCATE (WORK)
    1265           0 :                DEALLOCATE (data_copy)
    1266           0 :                DEALLOCATE (eigenvalues)
    1267             : 
    1268             :             END DO
    1269           0 :             CALL dbcsr_iterator_stop(iter)
    1270             : 
    1271           0 :             CALL dbcsr_finalize(matrix_t_blk_tmp)
    1272             :             CALL dbcsr_filter(matrix_t_blk_tmp, &
    1273           0 :                               almo_scf_env%eps_filter)
    1274             :             CALL dbcsr_copy(almo_scf_env%matrix_t_blk(ispin), &
    1275           0 :                             matrix_t_blk_tmp)
    1276           0 :             CALL dbcsr_release(matrix_t_blk_tmp)
    1277             : 
    1278             :          ELSE
    1279             : 
    1280             :             !! generate a random set of ALMOs
    1281             :             !! matrix_t_blk should already be initiated to the proper domain structure
    1282             :             CALL dbcsr_init_random(almo_scf_env%matrix_t_blk(ispin), &
    1283          60 :                                    keep_sparsity=.TRUE.)
    1284             : 
    1285             :             CALL dbcsr_create(matrix_t_blk_tmp, &
    1286             :                               template=almo_scf_env%matrix_t_blk(ispin), &
    1287          60 :                               matrix_type=dbcsr_type_no_symmetry)
    1288             : 
    1289             :             ! use current ALMOs in matrix_t_blk and project them onto the blocked dm
    1290             :             ! compute T_new = R_blk S_blk T_random
    1291             :             CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_s_blk(1), &
    1292             :                                 almo_scf_env%matrix_t_blk(ispin), &
    1293             :                                 0.0_dp, matrix_t_blk_tmp, &
    1294          60 :                                 filter_eps=almo_scf_env%eps_filter)
    1295             : 
    1296             :             CALL dbcsr_multiply("N", "N", 1.0_dp, &
    1297             :                                 almo_scf_env%matrix_p_blk(ispin), matrix_t_blk_tmp, &
    1298             :                                 0.0_dp, almo_scf_env%matrix_t_blk(ispin), &
    1299          60 :                                 filter_eps=almo_scf_env%eps_filter)
    1300             : 
    1301          60 :             CALL dbcsr_release(matrix_t_blk_tmp)
    1302             : 
    1303             :          END IF
    1304             : 
    1305             :       END DO
    1306             : 
    1307          60 :       CALL timestop(handle)
    1308             : 
    1309         120 :    END SUBROUTINE almo_scf_p_blk_to_t_blk
    1310             : 
    1311             : ! **************************************************************************************************
    1312             : !> \brief Apply an occupation-rescaling trick to ALMOs for smearing.
    1313             : !>        Partially occupied orbitals are considered full and rescaled by SQRT(occupation_number)
    1314             : !>        (this was designed to be used with smearing only)
    1315             : !> \param matrix_t ...
    1316             : !> \param mo_energies ...
    1317             : !> \param mu_of_domain ...
    1318             : !> \param real_ne_of_domain ...
    1319             : !> \param spin_kTS ...
    1320             : !> \param smear_e_temp ...
    1321             : !> \param ndomains ...
    1322             : !> \param nocc_of_domain ...
    1323             : !> \par History
    1324             : !>       2018.09 created [Ruben Staub]
    1325             : !> \author Ruben Staub
    1326             : ! **************************************************************************************************
    1327          40 :    SUBROUTINE almo_scf_t_rescaling(matrix_t, mo_energies, mu_of_domain, real_ne_of_domain, &
    1328          20 :                                    spin_kTS, smear_e_temp, ndomains, nocc_of_domain)
    1329             : 
    1330             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_t
    1331             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: mo_energies
    1332             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: mu_of_domain, real_ne_of_domain
    1333             :       REAL(KIND=dp), INTENT(INOUT)                       :: spin_kTS
    1334             :       REAL(KIND=dp), INTENT(IN)                          :: smear_e_temp
    1335             :       INTEGER, INTENT(IN)                                :: ndomains
    1336             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nocc_of_domain
    1337             : 
    1338             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'almo_scf_t_rescaling'
    1339             : 
    1340             :       INTEGER                                            :: handle, idomain, neigenval_used, nmo
    1341             :       REAL(KIND=dp)                                      :: kTS
    1342          20 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: occupation_numbers, rescaling_factors
    1343             : 
    1344          20 :       CALL timeset(routineN, handle)
    1345             : 
    1346             :       !!
    1347             :       !! Initialization
    1348             :       !!
    1349          20 :       nmo = SIZE(mo_energies)
    1350          60 :       ALLOCATE (occupation_numbers(nmo))
    1351          60 :       ALLOCATE (rescaling_factors(nmo))
    1352             : 
    1353             :       !!
    1354             :       !! Set occupation numbers for smearing
    1355             :       !!
    1356             :       !! occupation numbers are obtained using Fermi-Dirac smearing with orbital energies stored in mo_energies
    1357             :       !! nocc_of_domain is the number of partially occupied orbitals, while real_ne_of_domain is the real number of electrons
    1358          20 :       neigenval_used = 0 !! this is used as an offset to copy sub-arrays
    1359             : 
    1360             :       !! Reset electronic entropy term
    1361          20 :       spin_kTS = 0.0_dp
    1362             : 
    1363             :       !! Apply Fermi-Dirac smearing for each domain and store associated occupations for the whole system
    1364          60 :       DO idomain = 1, ndomains
    1365             :          CALL FermiFixed(occupation_numbers(1 + neigenval_used:nocc_of_domain(idomain) + neigenval_used), &
    1366             :                          mu_of_domain(idomain), &
    1367             :                          kTS, &
    1368             :                          mo_energies(1 + neigenval_used:nocc_of_domain(idomain) + neigenval_used), &
    1369             :                          real_ne_of_domain(idomain), &
    1370             :                          smear_e_temp, &
    1371          40 :                          1.0_dp) !! Warning, maxocc is set to 1 since we don't want to interfere with the spin_factor rescaling
    1372          40 :          spin_kTS = spin_kTS + kTS !! Add up electronic entropy contributions
    1373          60 :          neigenval_used = neigenval_used + nocc_of_domain(idomain) !! Update eigenvalues index offset
    1374             :       END DO
    1375         700 :       rescaling_factors(:) = SQRT(occupation_numbers) !! scale = sqrt(occupation_number)
    1376             : 
    1377             :       !!
    1378             :       !! Rescaling electronic entropy contribution by spin_factor (deprecated)
    1379             :       !! (currently, entropy is rescaled by spin_factor with the density matrix)
    1380             :       !!
    1381             :       !!IF (almo_scf_env%nspins == 1) THEN
    1382             :       !!   spin_kTS = spin_kTS*2.0_dp
    1383             :       !!END IF
    1384             : 
    1385             :       !!
    1386             :       !! Rescaling of T (i.e. ALMOs)
    1387             :       !!
    1388          20 :       CALL dbcsr_scale_by_vector(matrix_t, rescaling_factors, side='right') !! Apply occupation-rescaling trick
    1389             : 
    1390             :       !!
    1391             :       !! Debug tools (for debug purpose only)
    1392             :       !!
    1393             :       !! WRITE (*,*) "occupations", occupation_numbers(:) !! debug
    1394             :       !! WRITE (*,*) "eigenvalues", mo_energies(:) !! debug
    1395             :       !! WRITE (*,*) "kTS (spin_factor excluded) = ", spin_kTS !! debug
    1396             : 
    1397             :       !!
    1398             :       !! Cleaning up before exit
    1399             :       !!
    1400          20 :       DEALLOCATE (occupation_numbers)
    1401          20 :       DEALLOCATE (rescaling_factors)
    1402             : 
    1403          20 :       CALL timestop(handle)
    1404             : 
    1405          40 :    END SUBROUTINE almo_scf_t_rescaling
    1406             : 
    1407             : ! **************************************************************************************************
    1408             : !> \brief Computes the overlap matrix of MO orbitals
    1409             : !> \param bra ...
    1410             : !> \param ket ...
    1411             : !> \param overlap ...
    1412             : !> \param metric ...
    1413             : !> \param retain_overlap_sparsity ...
    1414             : !> \param eps_filter ...
    1415             : !> \param smear ...
    1416             : !> \par History
    1417             : !>       2011.08 created [Rustam Z Khaliullin]
    1418             : !>       2018.09 smearing support [Ruben Staub]
    1419             : !> \author Rustam Z Khaliullin
    1420             : ! **************************************************************************************************
    1421         378 :    SUBROUTINE get_overlap(bra, ket, overlap, metric, retain_overlap_sparsity, &
    1422             :                           eps_filter, smear)
    1423             : 
    1424             :       TYPE(dbcsr_type), INTENT(IN)                       :: bra, ket
    1425             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: overlap
    1426             :       TYPE(dbcsr_type), INTENT(IN)                       :: metric
    1427             :       LOGICAL, INTENT(IN), OPTIONAL                      :: retain_overlap_sparsity
    1428             :       REAL(KIND=dp)                                      :: eps_filter
    1429             :       LOGICAL, INTENT(IN), OPTIONAL                      :: smear
    1430             : 
    1431             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_overlap'
    1432             : 
    1433             :       INTEGER                                            :: dim0, handle
    1434             :       LOGICAL                                            :: local_retain_sparsity, smearing
    1435         378 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: diag_correction
    1436             :       TYPE(dbcsr_type)                                   :: tmp
    1437             : 
    1438         378 :       CALL timeset(routineN, handle)
    1439             : 
    1440         378 :       IF (.NOT. PRESENT(retain_overlap_sparsity)) THEN
    1441           0 :          local_retain_sparsity = .FALSE.
    1442             :       ELSE
    1443         378 :          local_retain_sparsity = retain_overlap_sparsity
    1444             :       END IF
    1445             : 
    1446         378 :       IF (.NOT. PRESENT(smear)) THEN
    1447             :          smearing = .FALSE.
    1448             :       ELSE
    1449         378 :          smearing = smear
    1450             :       END IF
    1451             : 
    1452             :       CALL dbcsr_create(tmp, template=ket, &
    1453         378 :                         matrix_type=dbcsr_type_no_symmetry)
    1454             : 
    1455             :       ! TMP=metric*ket
    1456             :       CALL dbcsr_multiply("N", "N", 1.0_dp, &
    1457             :                           metric, ket, 0.0_dp, tmp, &
    1458         378 :                           filter_eps=eps_filter)
    1459             : 
    1460             :       ! OVERLAP=tr(bra)*TMP
    1461             :       CALL dbcsr_multiply("T", "N", 1.0_dp, &
    1462             :                           bra, tmp, 0.0_dp, overlap, &
    1463             :                           retain_sparsity=local_retain_sparsity, &
    1464         378 :                           filter_eps=eps_filter)
    1465             : 
    1466         378 :       CALL dbcsr_release(tmp)
    1467             : 
    1468             :       !! If smearing ALMO is requested, apply correction of the occupation-rescaling trick
    1469             :       !! (i.e. converting rescaled orbitals into selfish orbitals)
    1470             :       !! (i.e. the diagonal blocks of the rescaled overlap must remain unscaled)
    1471             :       !! Since we have here orthonormal MOs within a fragment, diagonal blocks are identity matrices
    1472             :       !! Therefore, one only need to restore the diagonal to 1
    1473             :       !! RS-WARNING: Assume orthonormal MOs within a fragment
    1474         378 :       IF (smearing) THEN
    1475           4 :          CALL dbcsr_get_info(overlap, nfullrows_total=dim0)
    1476          12 :          ALLOCATE (diag_correction(dim0))
    1477         132 :          diag_correction = 1.0_dp
    1478           4 :          CALL dbcsr_set_diag(overlap, diag_correction)
    1479           8 :          DEALLOCATE (diag_correction)
    1480             :       END IF
    1481             : 
    1482         378 :       CALL timestop(handle)
    1483             : 
    1484         756 :    END SUBROUTINE get_overlap
    1485             : 
    1486             : ! **************************************************************************************************
    1487             : !> \brief orthogonalize MOs
    1488             : !> \param ket ...
    1489             : !> \param overlap ...
    1490             : !> \param metric ...
    1491             : !> \param retain_locality ...
    1492             : !> \param only_normalize ...
    1493             : !> \param nocc_of_domain ...
    1494             : !> \param eps_filter ...
    1495             : !> \param order_lanczos ...
    1496             : !> \param eps_lanczos ...
    1497             : !> \param max_iter_lanczos ...
    1498             : !> \param overlap_sqrti ...
    1499             : !> \param smear ...
    1500             : !> \par History
    1501             : !>       2012.03 created [Rustam Z Khaliullin]
    1502             : !>       2018.09 smearing support [Ruben Staub]
    1503             : !> \author Rustam Z Khaliullin
    1504             : ! **************************************************************************************************
    1505         378 :    SUBROUTINE orthogonalize_mos(ket, overlap, metric, retain_locality, only_normalize, &
    1506         378 :                                 nocc_of_domain, eps_filter, order_lanczos, eps_lanczos, &
    1507             :                                 max_iter_lanczos, overlap_sqrti, smear)
    1508             : 
    1509             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: ket, overlap
    1510             :       TYPE(dbcsr_type), INTENT(IN)                       :: metric
    1511             :       LOGICAL, INTENT(IN)                                :: retain_locality, only_normalize
    1512             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nocc_of_domain
    1513             :       REAL(KIND=dp)                                      :: eps_filter
    1514             :       INTEGER, INTENT(IN)                                :: order_lanczos
    1515             :       REAL(KIND=dp), INTENT(IN)                          :: eps_lanczos
    1516             :       INTEGER, INTENT(IN)                                :: max_iter_lanczos
    1517             :       TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL          :: overlap_sqrti
    1518             :       LOGICAL, INTENT(IN), OPTIONAL                      :: smear
    1519             : 
    1520             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'orthogonalize_mos'
    1521             : 
    1522             :       INTEGER                                            :: dim0, handle
    1523             :       LOGICAL                                            :: smearing
    1524         378 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: diagonal
    1525             :       TYPE(dbcsr_type)                                   :: matrix_sigma_blk_sqrt, &
    1526             :                                                             matrix_sigma_blk_sqrt_inv, &
    1527             :                                                             matrix_t_blk_tmp
    1528             : 
    1529         378 :       CALL timeset(routineN, handle)
    1530             : 
    1531         378 :       IF (.NOT. PRESENT(smear)) THEN
    1532         284 :          smearing = .FALSE.
    1533             :       ELSE
    1534          94 :          smearing = smear
    1535             :       END IF
    1536             : 
    1537             :       ! create block-diagonal sparsity pattern for the overlap
    1538             :       ! in case retain_locality is set to true
    1539             :       ! RZK-warning this will fail if distribution blocks are smaller than domains!!!
    1540         378 :       CALL dbcsr_set(overlap, 0.0_dp)
    1541         378 :       CALL dbcsr_add_on_diag(overlap, 1.0_dp)
    1542         378 :       CALL dbcsr_filter(overlap, eps_filter)
    1543             : 
    1544             :       CALL get_overlap(ket, ket, overlap, metric, retain_locality, &
    1545         378 :                        eps_filter, smear=smearing)
    1546             : 
    1547         378 :       IF (only_normalize) THEN
    1548             : 
    1549           0 :          CALL dbcsr_get_info(overlap, nfullrows_total=dim0)
    1550           0 :          ALLOCATE (diagonal(dim0))
    1551           0 :          CALL dbcsr_get_diag(overlap, diagonal)
    1552           0 :          CALL dbcsr_set(overlap, 0.0_dp)
    1553           0 :          CALL dbcsr_set_diag(overlap, diagonal)
    1554           0 :          DEALLOCATE (diagonal)
    1555           0 :          CALL dbcsr_filter(overlap, eps_filter)
    1556             : 
    1557             :       END IF
    1558             : 
    1559             :       CALL dbcsr_create(matrix_sigma_blk_sqrt, template=overlap, &
    1560         378 :                         matrix_type=dbcsr_type_no_symmetry)
    1561             :       CALL dbcsr_create(matrix_sigma_blk_sqrt_inv, template=overlap, &
    1562         378 :                         matrix_type=dbcsr_type_no_symmetry)
    1563             : 
    1564             :       ! compute sqrt and sqrt_inv of the blocked MO overlap
    1565         378 :       CALL set_zero_electron_blocks_in_mo_mo_matrix(overlap, nocc_of_domain, 1.0_dp)
    1566             :       CALL matrix_sqrt_Newton_Schulz(matrix_sigma_blk_sqrt, matrix_sigma_blk_sqrt_inv, &
    1567             :                                      overlap, threshold=eps_filter, &
    1568             :                                      order=order_lanczos, &
    1569             :                                      eps_lanczos=eps_lanczos, &
    1570         378 :                                      max_iter_lanczos=max_iter_lanczos)
    1571         378 :       CALL set_zero_electron_blocks_in_mo_mo_matrix(overlap, nocc_of_domain, 0.0_dp)
    1572             :       !CALL set_zero_electron_blocks_in_mo_mo_matrix(matrix_sigma_blk_sqrt,nocc_of_domain,0.0_dp)
    1573         378 :       CALL set_zero_electron_blocks_in_mo_mo_matrix(matrix_sigma_blk_sqrt_inv, nocc_of_domain, 0.0_dp)
    1574             : 
    1575             :       CALL dbcsr_create(matrix_t_blk_tmp, &
    1576             :                         template=ket, &
    1577         378 :                         matrix_type=dbcsr_type_no_symmetry)
    1578             : 
    1579             :       CALL dbcsr_multiply("N", "N", 1.0_dp, &
    1580             :                           ket, &
    1581             :                           matrix_sigma_blk_sqrt_inv, &
    1582             :                           0.0_dp, matrix_t_blk_tmp, &
    1583         378 :                           filter_eps=eps_filter)
    1584             : 
    1585             :       ! update the orbitals with the orthonormalized MOs
    1586         378 :       CALL dbcsr_copy(ket, matrix_t_blk_tmp)
    1587             : 
    1588             :       ! return overlap SQRT_INV if necessary
    1589         378 :       IF (PRESENT(overlap_sqrti)) THEN
    1590             :          CALL dbcsr_copy(overlap_sqrti, &
    1591           0 :                          matrix_sigma_blk_sqrt_inv)
    1592             :       END IF
    1593             : 
    1594         378 :       CALL dbcsr_release(matrix_t_blk_tmp)
    1595         378 :       CALL dbcsr_release(matrix_sigma_blk_sqrt)
    1596         378 :       CALL dbcsr_release(matrix_sigma_blk_sqrt_inv)
    1597             : 
    1598         378 :       CALL timestop(handle)
    1599             : 
    1600         756 :    END SUBROUTINE orthogonalize_mos
    1601             : 
    1602             : ! **************************************************************************************************
    1603             : !> \brief computes the idempotent density matrix from MOs
    1604             : !>        MOs can be either orthogonal or non-orthogonal
    1605             : !> \param t ...
    1606             : !> \param p ...
    1607             : !> \param eps_filter ...
    1608             : !> \param orthog_orbs ...
    1609             : !> \param nocc_of_domain ...
    1610             : !> \param s ...
    1611             : !> \param sigma ...
    1612             : !> \param sigma_inv ...
    1613             : !> \param use_guess ...
    1614             : !> \param smear ...
    1615             : !> \param algorithm to inver sigma: 0 - Hotelling (linear), 1 - Cholesky (cubic, low prefactor)
    1616             : !> \param para_env ...
    1617             : !> \param blacs_env ...
    1618             : !> \param eps_lanczos ...
    1619             : !> \param max_iter_lanczos ...
    1620             : !> \param inverse_accelerator ...
    1621             : !> \param inv_eps_factor ...
    1622             : !> \par History
    1623             : !>       2011.07 created [Rustam Z Khaliullin]
    1624             : !>       2018.09 smearing support [Ruben Staub]
    1625             : !> \author Rustam Z Khaliullin
    1626             : ! **************************************************************************************************
    1627        1948 :    SUBROUTINE almo_scf_t_to_proj(t, p, eps_filter, orthog_orbs, nocc_of_domain, s, sigma, sigma_inv, &
    1628             :                                  use_guess, smear, algorithm, para_env, blacs_env, eps_lanczos, &
    1629             :                                  max_iter_lanczos, inverse_accelerator, inv_eps_factor)
    1630             : 
    1631             :       TYPE(dbcsr_type), INTENT(IN)                       :: t
    1632             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: p
    1633             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
    1634             :       LOGICAL, INTENT(IN)                                :: orthog_orbs
    1635             :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: nocc_of_domain
    1636             :       TYPE(dbcsr_type), INTENT(IN), OPTIONAL             :: s
    1637             :       TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL          :: sigma, sigma_inv
    1638             :       LOGICAL, INTENT(IN), OPTIONAL                      :: use_guess, smear
    1639             :       INTEGER, INTENT(IN), OPTIONAL                      :: algorithm
    1640             :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env
    1641             :       TYPE(cp_blacs_env_type), OPTIONAL, POINTER         :: blacs_env
    1642             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: eps_lanczos
    1643             :       INTEGER, INTENT(IN), OPTIONAL                      :: max_iter_lanczos, inverse_accelerator
    1644             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: inv_eps_factor
    1645             : 
    1646             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'almo_scf_t_to_proj'
    1647             : 
    1648             :       INTEGER                                            :: dim0, handle, my_accelerator, &
    1649             :                                                             my_algorithm
    1650             :       LOGICAL                                            :: smearing, use_sigma_inv_guess
    1651             :       REAL(KIND=dp)                                      :: my_inv_eps_factor
    1652        1948 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: diag_correction
    1653             :       TYPE(dbcsr_type)                                   :: t_tmp
    1654             : 
    1655        1948 :       CALL timeset(routineN, handle)
    1656             : 
    1657             :       ! make sure that S, sigma and sigma_inv are present for non-orthogonal orbitals
    1658        1948 :       IF (.NOT. orthog_orbs) THEN
    1659             :          IF ((.NOT. PRESENT(s)) .OR. (.NOT. PRESENT(sigma)) .OR. &
    1660        1948 :              (.NOT. PRESENT(sigma_inv)) .OR. (.NOT. PRESENT(nocc_of_domain))) THEN
    1661           0 :             CPABORT("Nonorthogonal orbitals need more input")
    1662             :          END IF
    1663             :       END IF
    1664             : 
    1665        1948 :       my_algorithm = 0
    1666        1948 :       IF (PRESENT(algorithm)) my_algorithm = algorithm
    1667             : 
    1668        1948 :       IF (my_algorithm == 1 .AND. (.NOT. PRESENT(para_env) .OR. .NOT. PRESENT(blacs_env))) &
    1669           0 :          CPABORT("PARA and BLACS env are necessary for cholesky algorithm")
    1670             : 
    1671        1948 :       use_sigma_inv_guess = .FALSE.
    1672        1948 :       IF (PRESENT(use_guess)) THEN
    1673        1948 :          use_sigma_inv_guess = use_guess
    1674             :       END IF
    1675             : 
    1676        1948 :       IF (.NOT. PRESENT(smear)) THEN
    1677             :          smearing = .FALSE.
    1678             :       ELSE
    1679         466 :          smearing = smear
    1680             :       END IF
    1681             : 
    1682        1948 :       my_accelerator = 1
    1683        1948 :       IF (PRESENT(inverse_accelerator)) my_accelerator = inverse_accelerator
    1684             : 
    1685        1948 :       my_inv_eps_factor = 10.0_dp
    1686        1948 :       IF (PRESENT(inv_eps_factor)) my_inv_eps_factor = inv_eps_factor
    1687             : 
    1688        1948 :       IF (orthog_orbs) THEN
    1689             : 
    1690             :          CALL dbcsr_multiply("N", "T", 1.0_dp, t, t, &
    1691           0 :                              0.0_dp, p, filter_eps=eps_filter)
    1692             : 
    1693             :       ELSE
    1694             : 
    1695        1948 :          CALL dbcsr_create(t_tmp, template=t)
    1696             : 
    1697             :          ! TMP=S.T
    1698             :          CALL dbcsr_multiply("N", "N", 1.0_dp, s, t, 0.0_dp, t_tmp, &
    1699        1948 :                              filter_eps=eps_filter)
    1700             : 
    1701             :          ! Sig=tr(T).TMP - get MO overlap
    1702             :          CALL dbcsr_multiply("T", "N", 1.0_dp, t, t_tmp, 0.0_dp, sigma, &
    1703        1948 :                              filter_eps=eps_filter)
    1704             : 
    1705             :          !! If smearing ALMO is requested, apply correction of the occupation-rescaling trick
    1706             :          !! (i.e. converting rescaled orbitals into selfish orbitals)
    1707             :          !! (i.e. the diagonal blocks of the rescaled overlap must remain unscaled)
    1708             :          !! Since we have here orthonormal MOs within a fragment, diagonal blocks are identity matrices
    1709             :          !! Therefore, one only need to restore the diagonal to 1
    1710             :          !! RS-WARNING: Assume orthonormal MOs within a fragment
    1711        1948 :          IF (smearing) THEN
    1712          20 :             CALL dbcsr_get_info(sigma, nfullrows_total=dim0)
    1713          60 :             ALLOCATE (diag_correction(dim0))
    1714         700 :             diag_correction = 1.0_dp
    1715          20 :             CALL dbcsr_set_diag(sigma, diag_correction)
    1716          40 :             DEALLOCATE (diag_correction)
    1717             :          END IF
    1718             : 
    1719             :          ! invert MO overlap
    1720        1948 :          CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma, nocc_of_domain, 1.0_dp)
    1721          26 :          SELECT CASE (my_algorithm)
    1722             :          CASE (spd_inversion_ls_taylor)
    1723             : 
    1724             :             CALL invert_Taylor( &
    1725             :                matrix_inverse=sigma_inv, &
    1726             :                matrix=sigma, &
    1727             :                use_inv_as_guess=use_sigma_inv_guess, &
    1728             :                threshold=eps_filter*my_inv_eps_factor, &
    1729             :                filter_eps=eps_filter, &
    1730             :                !accelerator_order=my_accelerator, &
    1731             :                !eps_lanczos=eps_lanczos, &
    1732             :                !max_iter_lanczos=max_iter_lanczos, &
    1733          26 :                silent=.FALSE.)
    1734             : 
    1735             :          CASE (spd_inversion_ls_hotelling)
    1736             : 
    1737             :             CALL invert_Hotelling( &
    1738             :                matrix_inverse=sigma_inv, &
    1739             :                matrix=sigma, &
    1740             :                use_inv_as_guess=use_sigma_inv_guess, &
    1741             :                threshold=eps_filter*my_inv_eps_factor, &
    1742             :                filter_eps=eps_filter, &
    1743             :                accelerator_order=my_accelerator, &
    1744             :                eps_lanczos=eps_lanczos, &
    1745             :                max_iter_lanczos=max_iter_lanczos, &
    1746        1436 :                silent=.FALSE.)
    1747             : 
    1748             :          CASE (spd_inversion_dense_cholesky)
    1749             : 
    1750             :             ! invert using cholesky
    1751         486 :             CALL dbcsr_copy(sigma_inv, sigma)
    1752             :             CALL cp_dbcsr_cholesky_decompose(sigma_inv, &
    1753             :                                              para_env=para_env, &
    1754         486 :                                              blacs_env=blacs_env)
    1755             :             CALL cp_dbcsr_cholesky_invert(sigma_inv, &
    1756             :                                           para_env=para_env, &
    1757             :                                           blacs_env=blacs_env, &
    1758         486 :                                           upper_to_full=.TRUE.)
    1759         486 :             CALL dbcsr_filter(sigma_inv, eps_filter)
    1760             :          CASE DEFAULT
    1761        1948 :             CPABORT("Illegal MO overalp inversion algorithm")
    1762             :          END SELECT
    1763        1948 :          CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma, nocc_of_domain, 0.0_dp)
    1764        1948 :          CALL set_zero_electron_blocks_in_mo_mo_matrix(sigma_inv, nocc_of_domain, 0.0_dp)
    1765             : 
    1766             :          ! TMP=T.SigInv
    1767             :          CALL dbcsr_multiply("N", "N", 1.0_dp, t, sigma_inv, 0.0_dp, t_tmp, &
    1768        1948 :                              filter_eps=eps_filter)
    1769             : 
    1770             :          ! P=TMP.tr(T_blk)
    1771             :          CALL dbcsr_multiply("N", "T", 1.0_dp, t_tmp, t, 0.0_dp, p, &
    1772        1948 :                              filter_eps=eps_filter)
    1773             : 
    1774        1948 :          CALL dbcsr_release(t_tmp)
    1775             : 
    1776             :       END IF
    1777             : 
    1778        1948 :       CALL timestop(handle)
    1779             : 
    1780        3896 :    END SUBROUTINE almo_scf_t_to_proj
    1781             : 
    1782             : ! **************************************************************************************************
    1783             : !> \brief self-explanatory
    1784             : !> \param matrix ...
    1785             : !> \param nocc_of_domain ...
    1786             : !> \param value ...
    1787             : !> \param
    1788             : !> \par History
    1789             : !>       2016.12 created [Rustam Z Khaliullin]
    1790             : !> \author Rustam Z Khaliullin
    1791             : ! **************************************************************************************************
    1792        6978 :    SUBROUTINE set_zero_electron_blocks_in_mo_mo_matrix(matrix, nocc_of_domain, value)
    1793             : 
    1794             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
    1795             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nocc_of_domain
    1796             :       REAL(KIND=dp), INTENT(IN)                          :: value
    1797             : 
    1798             :       INTEGER                                            :: hold, iblock_col, iblock_row, mynode, &
    1799             :                                                             nblkrows_tot, row
    1800             :       LOGICAL                                            :: found, tr
    1801        6978 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_new_block
    1802             :       TYPE(dbcsr_distribution_type)                      :: dist
    1803             : 
    1804        6978 :       CALL dbcsr_get_info(matrix, distribution=dist)
    1805        6978 :       CALL dbcsr_distribution_get(dist, mynode=mynode)
    1806             :       !mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(dbcsr_distribution(matrix)))
    1807        6978 :       CALL dbcsr_work_create(matrix, work_mutable=.TRUE.)
    1808             : 
    1809        6978 :       nblkrows_tot = dbcsr_nblkrows_total(matrix)
    1810             : 
    1811       53406 :       DO row = 1, nblkrows_tot
    1812       53406 :          IF (nocc_of_domain(row) == 0) THEN
    1813        7632 :             tr = .FALSE.
    1814        7632 :             iblock_row = row
    1815        7632 :             iblock_col = row
    1816        7632 :             CALL dbcsr_get_stored_coordinates(matrix, iblock_row, iblock_col, hold)
    1817        7632 :             IF (hold .EQ. mynode) THEN
    1818        3816 :                NULLIFY (p_new_block)
    1819        3816 :                CALL dbcsr_get_block_p(matrix, iblock_row, iblock_col, p_new_block, found)
    1820        3816 :                IF (found) THEN
    1821        2792 :                   p_new_block(1, 1) = value
    1822             :                ELSE
    1823        1024 :                   CALL dbcsr_reserve_block2d(matrix, iblock_row, iblock_col, p_new_block)
    1824        1024 :                   CPASSERT(ASSOCIATED(p_new_block))
    1825        1024 :                   p_new_block(1, 1) = value
    1826             :                END IF
    1827             :             END IF ! mynode
    1828             :          END IF !zero-electron block
    1829             :       END DO
    1830             : 
    1831        6978 :       CALL dbcsr_finalize(matrix)
    1832             : 
    1833        6978 :    END SUBROUTINE set_zero_electron_blocks_in_mo_mo_matrix
    1834             : 
    1835             : ! **************************************************************************************************
    1836             : !> \brief applies projector to the orbitals
    1837             : !>        |psi_out> = P |psi_in>   OR   |psi_out> = (1-P) |psi_in>,
    1838             : !>        where P = |psi_proj> (<psi_proj|psi_roj>)^{-1} <psi_proj|
    1839             : !> \param psi_in ...
    1840             : !> \param psi_out ...
    1841             : !> \param psi_projector ...
    1842             : !> \param metric ...
    1843             : !> \param project_out ...
    1844             : !> \param psi_projector_orthogonal ...
    1845             : !> \param proj_in_template ...
    1846             : !> \param eps_filter ...
    1847             : !> \param sig_inv_projector ...
    1848             : !> \param sig_inv_template ...
    1849             : !> \par History
    1850             : !>       2011.10 created [Rustam Z Khaliullin]
    1851             : !> \author Rustam Z Khaliullin
    1852             : ! **************************************************************************************************
    1853           0 :    SUBROUTINE apply_projector(psi_in, psi_out, psi_projector, metric, project_out, &
    1854             :                               psi_projector_orthogonal, proj_in_template, eps_filter, sig_inv_projector, &
    1855             :                               sig_inv_template)
    1856             : 
    1857             :       TYPE(dbcsr_type), INTENT(IN)                       :: psi_in
    1858             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: psi_out
    1859             :       TYPE(dbcsr_type), INTENT(IN)                       :: psi_projector, metric
    1860             :       LOGICAL, INTENT(IN)                                :: project_out, psi_projector_orthogonal
    1861             :       TYPE(dbcsr_type), INTENT(IN)                       :: proj_in_template
    1862             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
    1863             :       TYPE(dbcsr_type), INTENT(IN), OPTIONAL             :: sig_inv_projector, sig_inv_template
    1864             : 
    1865             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'apply_projector'
    1866             : 
    1867             :       INTEGER                                            :: handle
    1868             :       TYPE(dbcsr_type)                                   :: tmp_no, tmp_ov, tmp_ov2, tmp_sig, &
    1869             :                                                             tmp_sig_inv
    1870             : 
    1871           0 :       CALL timeset(routineN, handle)
    1872             : 
    1873             :       ! =S*PSI_proj
    1874           0 :       CALL dbcsr_create(tmp_no, template=psi_projector)
    1875             :       CALL dbcsr_multiply("N", "N", 1.0_dp, &
    1876             :                           metric, psi_projector, &
    1877             :                           0.0_dp, tmp_no, &
    1878           0 :                           filter_eps=eps_filter)
    1879             : 
    1880             :       ! =tr(S.PSI_proj)*PSI_in
    1881           0 :       CALL dbcsr_create(tmp_ov, template=proj_in_template)
    1882             :       CALL dbcsr_multiply("T", "N", 1.0_dp, &
    1883             :                           tmp_no, psi_in, &
    1884             :                           0.0_dp, tmp_ov, &
    1885           0 :                           filter_eps=eps_filter)
    1886             : 
    1887           0 :       IF (.NOT. psi_projector_orthogonal) THEN
    1888             :          ! =SigInv_proj*Sigma_OV
    1889             :          CALL dbcsr_create(tmp_ov2, &
    1890           0 :                            template=proj_in_template)
    1891           0 :          IF (PRESENT(sig_inv_projector)) THEN
    1892             :             CALL dbcsr_create(tmp_sig_inv, &
    1893           0 :                               template=sig_inv_projector)
    1894           0 :             CALL dbcsr_copy(tmp_sig_inv, sig_inv_projector)
    1895             :          ELSE
    1896           0 :             IF (.NOT. PRESENT(sig_inv_template)) THEN
    1897           0 :                CPABORT("PROGRAMMING ERROR: provide either template or sig_inv")
    1898             :             END IF
    1899             :             ! compute inverse overlap of the projector orbitals
    1900             :             CALL dbcsr_create(tmp_sig, &
    1901             :                               template=sig_inv_template, &
    1902           0 :                               matrix_type=dbcsr_type_no_symmetry)
    1903             :             CALL dbcsr_multiply("T", "N", 1.0_dp, &
    1904             :                                 psi_projector, tmp_no, 0.0_dp, tmp_sig, &
    1905           0 :                                 filter_eps=eps_filter)
    1906             :             CALL dbcsr_create(tmp_sig_inv, &
    1907             :                               template=sig_inv_template, &
    1908           0 :                               matrix_type=dbcsr_type_no_symmetry)
    1909             :             CALL invert_Hotelling(tmp_sig_inv, tmp_sig, &
    1910           0 :                                   threshold=eps_filter)
    1911           0 :             CALL dbcsr_release(tmp_sig)
    1912             :          END IF
    1913             :          CALL dbcsr_multiply("N", "N", 1.0_dp, &
    1914             :                              tmp_sig_inv, tmp_ov, 0.0_dp, tmp_ov2, &
    1915           0 :                              filter_eps=eps_filter)
    1916           0 :          CALL dbcsr_release(tmp_sig_inv)
    1917           0 :          CALL dbcsr_copy(tmp_ov, tmp_ov2)
    1918           0 :          CALL dbcsr_release(tmp_ov2)
    1919             :       END IF
    1920           0 :       CALL dbcsr_release(tmp_no)
    1921             : 
    1922             :       ! =PSI_proj*TMP_OV
    1923             :       CALL dbcsr_multiply("N", "N", 1.0_dp, &
    1924             :                           psi_projector, tmp_ov, 0.0_dp, psi_out, &
    1925           0 :                           filter_eps=eps_filter)
    1926           0 :       CALL dbcsr_release(tmp_ov)
    1927             : 
    1928             :       ! V_out=V_in-V_out
    1929           0 :       IF (project_out) THEN
    1930           0 :          CALL dbcsr_add(psi_out, psi_in, -1.0_dp, +1.0_dp)
    1931             :       END IF
    1932             : 
    1933           0 :       CALL timestop(handle)
    1934             : 
    1935           0 :    END SUBROUTINE apply_projector
    1936             : 
    1937             : !! **************************************************************************************************
    1938             : !!> \brief projects the occupied space out from the provided orbitals
    1939             : !!> \par History
    1940             : !!>       2011.07 created [Rustam Z Khaliullin]
    1941             : !!> \author Rustam Z Khaliullin
    1942             : !! **************************************************************************************************
    1943             : !  SUBROUTINE almo_scf_p_out_from_v(v_in,v_out,ov_template,ispin,almo_scf_env)
    1944             : !
    1945             : !    TYPE(dbcsr_type), INTENT(IN)                :: v_in, ov_template
    1946             : !    TYPE(dbcsr_type), INTENT(INOUT)             :: v_out
    1947             : !    INTEGER, INTENT(IN)                            :: ispin
    1948             : !    TYPE(almo_scf_env_type), INTENT(INOUT)         :: almo_scf_env
    1949             : !
    1950             : !    CHARACTER(LEN=*), PARAMETER :: &
    1951             : !      routineN = 'almo_scf_p_out_from_v', &
    1952             : !      routineP = moduleN//':'//routineN
    1953             : !
    1954             : !    TYPE(dbcsr_type)                      :: tmp_on, tmp_ov, tmp_ov2
    1955             : !    INTEGER                                  :: handle
    1956             : !    LOGICAL                                  :: failure
    1957             : !
    1958             : !    CALL timeset(routineN,handle)
    1959             : !
    1960             : !    ! =tr(T_blk)*S
    1961             : !    CALL dbcsr_init(tmp_on)
    1962             : !    CALL dbcsr_create(tmp_on,&
    1963             : !            template=almo_scf_env%matrix_t_tr(ispin))
    1964             : !    CALL dbcsr_multiply("T","N",1.0_dp,&
    1965             : !            almo_scf_env%matrix_t_blk(ispin),&
    1966             : !            almo_scf_env%matrix_s(1),&
    1967             : !            0.0_dp,tmp_on,&
    1968             : !            filter_eps=almo_scf_env%eps_filter)
    1969             : !
    1970             : !    ! =tr(T_blk).S*V_in
    1971             : !    CALL dbcsr_init(tmp_ov)
    1972             : !    CALL dbcsr_create(tmp_ov,template=ov_template)
    1973             : !    CALL dbcsr_multiply("N","N",1.0_dp,&
    1974             : !            tmp_on,v_in,0.0_dp,tmp_ov,&
    1975             : !            filter_eps=almo_scf_env%eps_filter)
    1976             : !    CALL dbcsr_release(tmp_on)
    1977             : !
    1978             : !    ! =SigmaInv*Sigma_OV
    1979             : !    CALL dbcsr_init(tmp_ov2)
    1980             : !    CALL dbcsr_create(tmp_ov2,template=ov_template)
    1981             : !    CALL dbcsr_multiply("N","N",1.0_dp,&
    1982             : !            almo_scf_env%matrix_sigma_inv(ispin),&
    1983             : !            tmp_ov,0.0_dp,tmp_ov2,&
    1984             : !            filter_eps=almo_scf_env%eps_filter)
    1985             : !    CALL dbcsr_release(tmp_ov)
    1986             : !
    1987             : !    ! =T_blk*SigmaInv.Sigma_OV
    1988             : !    CALL dbcsr_multiply("N","N",1.0_dp,&
    1989             : !            almo_scf_env%matrix_t_blk(ispin),&
    1990             : !            tmp_ov2,0.0_dp,v_out,&
    1991             : !            filter_eps=almo_scf_env%eps_filter)
    1992             : !    CALL dbcsr_release(tmp_ov2)
    1993             : !
    1994             : !    ! V_out=V_in-V_out=
    1995             : !    CALL dbcsr_add(v_out,v_in,-1.0_dp,+1.0_dp)
    1996             : !
    1997             : !    CALL timestop(handle)
    1998             : !
    1999             : !  END SUBROUTINE almo_scf_p_out_from_v
    2000             : 
    2001             : ! **************************************************************************************************
    2002             : !> \brief computes a unitary matrix from an arbitrary "generator" matrix
    2003             : !>        U = ( 1 - X + tr(X) ) ( 1 + X - tr(X) )^(-1)
    2004             : !> \param X ...
    2005             : !> \param U ...
    2006             : !> \param eps_filter ...
    2007             : !> \par History
    2008             : !>       2011.08 created [Rustam Z Khaliullin]
    2009             : !> \author Rustam Z Khaliullin
    2010             : ! **************************************************************************************************
    2011           0 :    SUBROUTINE generator_to_unitary(X, U, eps_filter)
    2012             : 
    2013             :       TYPE(dbcsr_type), INTENT(IN)                       :: X
    2014             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: U
    2015             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
    2016             : 
    2017             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'generator_to_unitary'
    2018             : 
    2019             :       INTEGER                                            :: handle, unit_nr
    2020             :       LOGICAL                                            :: safe_mode
    2021             :       REAL(KIND=dp)                                      :: frob_matrix, frob_matrix_base
    2022             :       TYPE(cp_logger_type), POINTER                      :: logger
    2023             :       TYPE(dbcsr_type)                                   :: delta, t1, t2, tmp1
    2024             : 
    2025           0 :       CALL timeset(routineN, handle)
    2026             : 
    2027           0 :       safe_mode = .TRUE.
    2028             : 
    2029             :       ! get a useful output_unit
    2030           0 :       logger => cp_get_default_logger()
    2031           0 :       IF (logger%para_env%is_source()) THEN
    2032           0 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    2033             :       ELSE
    2034             :          unit_nr = -1
    2035             :       END IF
    2036             : 
    2037             :       CALL dbcsr_create(t1, template=X, &
    2038           0 :                         matrix_type=dbcsr_type_no_symmetry)
    2039             :       CALL dbcsr_create(t2, template=X, &
    2040           0 :                         matrix_type=dbcsr_type_no_symmetry)
    2041             : 
    2042             :       ! create antisymmetric Delta = -X + tr(X)
    2043             :       CALL dbcsr_create(delta, template=X, &
    2044           0 :                         matrix_type=dbcsr_type_no_symmetry)
    2045           0 :       CALL dbcsr_transposed(delta, X)
    2046             : ! check that transposed is added correctly
    2047           0 :       CALL dbcsr_add(delta, X, 1.0_dp, -1.0_dp)
    2048             : 
    2049             :       ! compute (1 - Delta)^(-1)
    2050           0 :       CALL dbcsr_add_on_diag(t1, 1.0_dp)
    2051           0 :       CALL dbcsr_add(t1, delta, 1.0_dp, -1.0_dp)
    2052           0 :       CALL invert_Hotelling(t2, t1, threshold=eps_filter)
    2053             : 
    2054             :       IF (safe_mode) THEN
    2055             : 
    2056             :          CALL dbcsr_create(tmp1, template=X, &
    2057           0 :                            matrix_type=dbcsr_type_no_symmetry)
    2058             :          CALL dbcsr_multiply("N", "N", 1.0_dp, t2, t1, 0.0_dp, tmp1, &
    2059           0 :                              filter_eps=eps_filter)
    2060           0 :          frob_matrix_base = dbcsr_frobenius_norm(tmp1)
    2061           0 :          CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
    2062           0 :          frob_matrix = dbcsr_frobenius_norm(tmp1)
    2063           0 :          IF (unit_nr > 0) WRITE (unit_nr, *) "Error for (inv(A)*A-I)", frob_matrix/frob_matrix_base
    2064           0 :          CALL dbcsr_release(tmp1)
    2065             :       END IF
    2066             : 
    2067             :       CALL dbcsr_multiply("N", "N", 1.0_dp, delta, t2, 0.0_dp, U, &
    2068           0 :                           filter_eps=eps_filter)
    2069           0 :       CALL dbcsr_add(U, t2, 1.0_dp, 1.0_dp)
    2070             : 
    2071             :       IF (safe_mode) THEN
    2072             : 
    2073             :          CALL dbcsr_create(tmp1, template=X, &
    2074           0 :                            matrix_type=dbcsr_type_no_symmetry)
    2075             :          CALL dbcsr_multiply("T", "N", 1.0_dp, U, U, 0.0_dp, tmp1, &
    2076           0 :                              filter_eps=eps_filter)
    2077           0 :          frob_matrix_base = dbcsr_frobenius_norm(tmp1)
    2078           0 :          CALL dbcsr_add_on_diag(tmp1, -1.0_dp)
    2079           0 :          frob_matrix = dbcsr_frobenius_norm(tmp1)
    2080           0 :          IF (unit_nr > 0) WRITE (unit_nr, *) "Error for (trn(U)*U-I)", frob_matrix/frob_matrix_base
    2081           0 :          CALL dbcsr_release(tmp1)
    2082             :       END IF
    2083             : 
    2084           0 :       CALL timestop(handle)
    2085             : 
    2086           0 :    END SUBROUTINE generator_to_unitary
    2087             : 
    2088             : ! **************************************************************************************************
    2089             : !> \brief Parallel code for domain specific operations (my_action)
    2090             : !>         0. out = op1 * in
    2091             : !>         1. out = in - op2 * op1 * in
    2092             : !> \param matrix_in ...
    2093             : !> \param matrix_out ...
    2094             : !> \param operator1 ...
    2095             : !> \param operator2 ...
    2096             : !> \param dpattern ...
    2097             : !> \param map ...
    2098             : !> \param node_of_domain ...
    2099             : !> \param my_action ...
    2100             : !> \param filter_eps ...
    2101             : !> \param matrix_trimmer ...
    2102             : !> \param use_trimmer ...
    2103             : !> \par History
    2104             : !>       2013.01 created [Rustam Z. Khaliullin]
    2105             : !> \author Rustam Z. Khaliullin
    2106             : ! **************************************************************************************************
    2107        2394 :    SUBROUTINE apply_domain_operators(matrix_in, matrix_out, operator1, operator2, &
    2108        2394 :                                      dpattern, map, node_of_domain, my_action, filter_eps, matrix_trimmer, use_trimmer)
    2109             : 
    2110             :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_in
    2111             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_out
    2112             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2113             :          INTENT(IN)                                      :: operator1
    2114             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2115             :          INTENT(IN), OPTIONAL                            :: operator2
    2116             :       TYPE(dbcsr_type), INTENT(IN)                       :: dpattern
    2117             :       TYPE(domain_map_type), INTENT(IN)                  :: map
    2118             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: node_of_domain
    2119             :       INTEGER, INTENT(IN)                                :: my_action
    2120             :       REAL(KIND=dp)                                      :: filter_eps
    2121             :       TYPE(dbcsr_type), INTENT(IN), OPTIONAL             :: matrix_trimmer
    2122             :       LOGICAL, INTENT(IN), OPTIONAL                      :: use_trimmer
    2123             : 
    2124             :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_domain_operators'
    2125             : 
    2126             :       INTEGER                                            :: handle, ndomains
    2127             :       LOGICAL                                            :: matrix_trimmer_required, my_use_trimmer, &
    2128             :                                                             operator2_required
    2129             :       TYPE(domain_submatrix_type), ALLOCATABLE, &
    2130        2394 :          DIMENSION(:)                                    :: subm_in, subm_out, subm_temp
    2131             : 
    2132        2394 :       CALL timeset(routineN, handle)
    2133             : 
    2134        2394 :       my_use_trimmer = .FALSE.
    2135        2394 :       IF (PRESENT(use_trimmer)) THEN
    2136         612 :          my_use_trimmer = use_trimmer
    2137             :       END IF
    2138             : 
    2139        2394 :       operator2_required = .FALSE.
    2140        2394 :       matrix_trimmer_required = .FALSE.
    2141             : 
    2142        2394 :       IF (my_action .EQ. 1) operator2_required = .TRUE.
    2143             : 
    2144        2394 :       IF (my_use_trimmer) THEN
    2145           0 :          matrix_trimmer_required = .TRUE.
    2146           0 :          CPABORT("TRIMMED PROJECTOR DISABLED!")
    2147             :       END IF
    2148             : 
    2149        2394 :       IF (.NOT. PRESENT(operator2) .AND. operator2_required) THEN
    2150           0 :          CPABORT("SECOND OPERATOR IS REQUIRED")
    2151             :       END IF
    2152        2394 :       IF (.NOT. PRESENT(matrix_trimmer) .AND. matrix_trimmer_required) THEN
    2153           0 :          CPABORT("TRIMMER MATRIX IS REQUIRED")
    2154             :       END IF
    2155             : 
    2156        2394 :       ndomains = dbcsr_nblkcols_total(dpattern)
    2157             : 
    2158       22832 :       ALLOCATE (subm_in(ndomains))
    2159       22832 :       ALLOCATE (subm_temp(ndomains))
    2160       22832 :       ALLOCATE (subm_out(ndomains))
    2161             :       !!!TRIM ALLOCATE(subm_trimmer(ndomains))
    2162        2394 :       CALL init_submatrices(subm_in)
    2163        2394 :       CALL init_submatrices(subm_temp)
    2164        2394 :       CALL init_submatrices(subm_out)
    2165             : 
    2166             :       CALL construct_submatrices(matrix_in, subm_in, &
    2167        2394 :                                  dpattern, map, node_of_domain, select_row)
    2168             : 
    2169             :       !!!TRIM IF (matrix_trimmer_required) THEN
    2170             :       !!!TRIM    CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
    2171             :       !!!TRIM            dpattern,map,node_of_domain,select_row)
    2172             :       !!!TRIM ENDIF
    2173             : 
    2174        2394 :       IF (my_action .EQ. 0) THEN
    2175             :          ! for example, apply preconditioner
    2176             :          CALL multiply_submatrices('N', 'N', 1.0_dp, operator1, &
    2177        1552 :                                    subm_in, 0.0_dp, subm_out)
    2178         842 :       ELSE IF (my_action .EQ. 1) THEN
    2179             :          ! use for projectors
    2180         842 :          CALL copy_submatrices(subm_in, subm_out, .TRUE.)
    2181             :          CALL multiply_submatrices('N', 'N', 1.0_dp, operator1, &
    2182         842 :                                    subm_in, 0.0_dp, subm_temp)
    2183             :          CALL multiply_submatrices('N', 'N', -1.0_dp, operator2, &
    2184         842 :                                    subm_temp, 1.0_dp, subm_out)
    2185             :       ELSE
    2186           0 :          CPABORT("ILLEGAL ACTION")
    2187             :       END IF
    2188             : 
    2189        2394 :       CALL construct_dbcsr_from_submatrices(matrix_out, subm_out, dpattern)
    2190        2394 :       CALL dbcsr_filter(matrix_out, filter_eps)
    2191             : 
    2192        2394 :       CALL release_submatrices(subm_out)
    2193        2394 :       CALL release_submatrices(subm_temp)
    2194        2394 :       CALL release_submatrices(subm_in)
    2195             : 
    2196       18044 :       DEALLOCATE (subm_out)
    2197       18044 :       DEALLOCATE (subm_temp)
    2198       18044 :       DEALLOCATE (subm_in)
    2199             : 
    2200        2394 :       CALL timestop(handle)
    2201             : 
    2202        2394 :    END SUBROUTINE apply_domain_operators
    2203             : 
    2204             : ! **************************************************************************************************
    2205             : !> \brief Constructs preconditioners for each domain
    2206             : !>        -1. projected preconditioner
    2207             : !>         0. simple preconditioner
    2208             : !> \param matrix_main ...
    2209             : !> \param subm_s_inv ...
    2210             : !> \param subm_s_inv_half ...
    2211             : !> \param subm_s_half ...
    2212             : !> \param subm_r_down ...
    2213             : !> \param matrix_trimmer ...
    2214             : !> \param dpattern ...
    2215             : !> \param map ...
    2216             : !> \param node_of_domain ...
    2217             : !> \param preconditioner ...
    2218             : !> \param bad_modes_projector_down ...
    2219             : !> \param use_trimmer ...
    2220             : !> \param eps_zero_eigenvalues ...
    2221             : !> \param my_action ...
    2222             : !> \param skip_inversion ...
    2223             : !> \par History
    2224             : !>       2013.01 created [Rustam Z. Khaliullin]
    2225             : !> \author Rustam Z. Khaliullin
    2226             : ! **************************************************************************************************
    2227         550 :    SUBROUTINE construct_domain_preconditioner(matrix_main, subm_s_inv, &
    2228         550 :                                               subm_s_inv_half, subm_s_half, &
    2229         550 :                                               subm_r_down, matrix_trimmer, &
    2230         550 :                                               dpattern, map, node_of_domain, &
    2231         550 :                                               preconditioner, &
    2232         550 :                                               bad_modes_projector_down, &
    2233             :                                               use_trimmer, &
    2234             :                                               eps_zero_eigenvalues, &
    2235             :                                               my_action, skip_inversion)
    2236             : 
    2237             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_main
    2238             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2239             :          INTENT(IN), OPTIONAL                            :: subm_s_inv, subm_s_inv_half, &
    2240             :                                                             subm_s_half, subm_r_down
    2241             :       TYPE(dbcsr_type), INTENT(IN), OPTIONAL             :: matrix_trimmer
    2242             :       TYPE(dbcsr_type), INTENT(IN)                       :: dpattern
    2243             :       TYPE(domain_map_type), INTENT(IN)                  :: map
    2244             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: node_of_domain
    2245             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2246             :          INTENT(INOUT)                                   :: preconditioner
    2247             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2248             :          INTENT(INOUT), OPTIONAL                         :: bad_modes_projector_down
    2249             :       LOGICAL, INTENT(IN), OPTIONAL                      :: use_trimmer
    2250             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: eps_zero_eigenvalues
    2251             :       INTEGER, INTENT(IN)                                :: my_action
    2252             :       LOGICAL, INTENT(IN), OPTIONAL                      :: skip_inversion
    2253             : 
    2254             :       CHARACTER(len=*), PARAMETER :: routineN = 'construct_domain_preconditioner'
    2255             : 
    2256             :       INTEGER                                            :: handle, idomain, index1_end, &
    2257             :                                                             index1_start, n_domain_mos, naos, &
    2258             :                                                             nblkrows_tot, ndomains, neighbor, row
    2259         550 :       INTEGER, DIMENSION(:), POINTER                     :: nmos
    2260             :       LOGICAL :: eps_zero_eigenvalues_required, matrix_r_required, matrix_s_half_required, &
    2261             :          matrix_s_inv_half_required, matrix_s_inv_required, matrix_trimmer_required, &
    2262             :          my_skip_inversion, my_use_trimmer
    2263         550 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Minv, proj_array
    2264             :       TYPE(domain_submatrix_type), ALLOCATABLE, &
    2265         550 :          DIMENSION(:)                                    :: subm_main, subm_tmp, subm_tmp2
    2266             : 
    2267         550 :       CALL timeset(routineN, handle)
    2268             : 
    2269         550 :       my_use_trimmer = .FALSE.
    2270         550 :       IF (PRESENT(use_trimmer)) THEN
    2271         550 :          my_use_trimmer = use_trimmer
    2272             :       END IF
    2273             : 
    2274         550 :       my_skip_inversion = .FALSE.
    2275         550 :       IF (PRESENT(skip_inversion)) THEN
    2276         550 :          my_skip_inversion = skip_inversion
    2277             :       END IF
    2278             : 
    2279         550 :       matrix_s_inv_half_required = .FALSE.
    2280         550 :       matrix_s_half_required = .FALSE.
    2281         550 :       eps_zero_eigenvalues_required = .FALSE.
    2282         550 :       matrix_s_inv_required = .FALSE.
    2283         550 :       matrix_trimmer_required = .FALSE.
    2284         550 :       matrix_r_required = .FALSE.
    2285             : 
    2286         550 :       IF (my_action .EQ. -1) matrix_s_inv_required = .TRUE.
    2287             :       IF (my_action .EQ. -1) matrix_r_required = .TRUE.
    2288         550 :       IF (my_use_trimmer) THEN
    2289           0 :          matrix_trimmer_required = .TRUE.
    2290           0 :          CPABORT("TRIMMED PRECONDITIONER DISABLED!")
    2291             :       END IF
    2292             :       ! tie the following optional arguments together to prevent bad calls
    2293         550 :       IF (PRESENT(bad_modes_projector_down)) THEN
    2294          18 :          matrix_s_inv_half_required = .TRUE.
    2295          18 :          matrix_s_half_required = .TRUE.
    2296          18 :          eps_zero_eigenvalues_required = .TRUE.
    2297             :       END IF
    2298             : 
    2299             :       ! check if all required optional arguments are provided
    2300         550 :       IF (.NOT. PRESENT(subm_s_inv_half) .AND. matrix_s_inv_half_required) THEN
    2301           0 :          CPABORT("S_inv_half SUBMATRICES ARE REQUIRED")
    2302             :       END IF
    2303         550 :       IF (.NOT. PRESENT(subm_s_half) .AND. matrix_s_half_required) THEN
    2304           0 :          CPABORT("S_half SUBMATRICES ARE REQUIRED")
    2305             :       END IF
    2306         550 :       IF (.NOT. PRESENT(eps_zero_eigenvalues) .AND. eps_zero_eigenvalues_required) THEN
    2307           0 :          CPABORT("EPS_ZERO_EIGENVALUES IS REQUIRED")
    2308             :       END IF
    2309         550 :       IF (.NOT. PRESENT(subm_s_inv) .AND. matrix_s_inv_required) THEN
    2310           0 :          CPABORT("S_inv SUBMATRICES ARE REQUIRED")
    2311             :       END IF
    2312         550 :       IF (.NOT. PRESENT(subm_r_down) .AND. matrix_r_required) THEN
    2313           0 :          CPABORT("R SUBMATRICES ARE REQUIRED")
    2314             :       END IF
    2315         550 :       IF (.NOT. PRESENT(matrix_trimmer) .AND. matrix_trimmer_required) THEN
    2316           0 :          CPABORT("TRIMMER MATRIX IS REQUIRED")
    2317             :       END IF
    2318             : 
    2319             :       CALL dbcsr_get_info(dpattern, &
    2320             :                           nblkcols_total=ndomains, &
    2321             :                           nblkrows_total=nblkrows_tot, &
    2322         550 :                           col_blk_size=nmos)
    2323             : 
    2324        4516 :       ALLOCATE (subm_main(ndomains))
    2325         550 :       CALL init_submatrices(subm_main)
    2326             :       !!!TRIM ALLOCATE(subm_trimmer(ndomains))
    2327             : 
    2328             :       CALL construct_submatrices(matrix_main, subm_main, &
    2329         550 :                                  dpattern, map, node_of_domain, select_row_col)
    2330             : 
    2331             :       !!!TRIM IF (matrix_trimmer_required) THEN
    2332             :       !!!TRIM    CALL construct_submatrices(matrix_trimmer,subm_trimmer,&
    2333             :       !!!TRIM            dpattern,map,node_of_domain,select_row)
    2334             :       !!!TRIM ENDIF
    2335             : 
    2336         550 :       IF (my_action .EQ. -1) THEN
    2337             :          ! project out the local occupied space
    2338             :          !tmp=MATMUL(subm_r(idomain)%mdata,Minv)
    2339             :          !Minv=MATMUL(tmp,subm_main(idomain)%mdata)
    2340             :          !subm_main(idomain)%mdata=subm_main(idomain)%mdata-&
    2341             :          !   Minv-TRANSPOSE(Minv)+MATMUL(Minv,TRANSPOSE(tmp))
    2342         248 :          ALLOCATE (subm_tmp(ndomains))
    2343         248 :          ALLOCATE (subm_tmp2(ndomains))
    2344          26 :          CALL init_submatrices(subm_tmp)
    2345          26 :          CALL init_submatrices(subm_tmp2)
    2346             :          CALL multiply_submatrices('N', 'N', 1.0_dp, subm_r_down, &
    2347          26 :                                    subm_s_inv, 0.0_dp, subm_tmp)
    2348             :          CALL multiply_submatrices('N', 'N', 1.0_dp, subm_tmp, &
    2349          26 :                                    subm_main, 0.0_dp, subm_tmp2)
    2350          26 :          CALL add_submatrices(1.0_dp, subm_main, -1.0_dp, subm_tmp2, 'N')
    2351          26 :          CALL add_submatrices(1.0_dp, subm_main, -1.0_dp, subm_tmp2, 'T')
    2352             :          CALL multiply_submatrices('N', 'T', 1.0_dp, subm_tmp2, &
    2353          26 :                                    subm_tmp, 1.0_dp, subm_main)
    2354          26 :          CALL release_submatrices(subm_tmp)
    2355          26 :          CALL release_submatrices(subm_tmp2)
    2356         196 :          DEALLOCATE (subm_tmp2)
    2357         196 :          DEALLOCATE (subm_tmp)
    2358             :       END IF
    2359             : 
    2360         550 :       IF (my_skip_inversion) THEN
    2361         316 :          CALL copy_submatrices(subm_main, preconditioner, .TRUE.)
    2362             :       ELSE
    2363             :          ! loop over domains - perform inversion
    2364        1284 :          DO idomain = 1, ndomains
    2365             : 
    2366             :             ! check if the submatrix exists
    2367        1284 :             IF (subm_main(idomain)%domain .GT. 0) THEN
    2368             : 
    2369             :                ! find sizes of MO submatrices
    2370         525 :                IF (idomain .EQ. 1) THEN
    2371             :                   index1_start = 1
    2372             :                ELSE
    2373         408 :                   index1_start = map%index1(idomain - 1)
    2374             :                END IF
    2375         525 :                index1_end = map%index1(idomain) - 1
    2376             : 
    2377         525 :                n_domain_mos = 0
    2378        2320 :                DO row = index1_start, index1_end
    2379        1795 :                   neighbor = map%pairs(row, 1)
    2380        2320 :                   n_domain_mos = n_domain_mos + nmos(neighbor)
    2381             :                END DO
    2382             : 
    2383         525 :                naos = subm_main(idomain)%nrows
    2384             :                !WRITE(*,*) "Domain, mo_self_and_neig, ao_domain: ", idomain, n_domain_mos, naos
    2385             : 
    2386        2100 :                ALLOCATE (Minv(naos, naos))
    2387             : 
    2388             :                !!!TRIM IF (my_use_trimmer) THEN
    2389             :                !!!TRIM    ! THIS IS SUPER EXPENSIVE (ELIMINATE)
    2390             :                !!!TRIM    ! trim the main matrix before inverting
    2391             :                !!!TRIM    ! assume that the trimmer columns are different (i.e. the main matrix is different for each MO)
    2392             :                !!!TRIM    allocate(tmp(naos,nmos(idomain)))
    2393             :                !!!TRIM    DO ii=1, nmos(idomain)
    2394             :                !!!TRIM       ! transform the main matrix using the trimmer for the current MO
    2395             :                !!!TRIM       DO jj=1, naos
    2396             :                !!!TRIM          DO kk=1, naos
    2397             :                !!!TRIM             Mstore(jj,kk)=sumb_main(idomain)%mdata(jj,kk)*&
    2398             :                !!!TRIM                subm_trimmer(idomain)%mdata(jj,ii)*&
    2399             :                !!!TRIM                subm_trimmer(idomain)%mdata(kk,ii)
    2400             :                !!!TRIM          ENDDO
    2401             :                !!!TRIM       ENDDO
    2402             :                !!!TRIM       ! invert the main matrix (exclude some eigenvalues, shift some)
    2403             :                !!!TRIM       CALL pseudo_invert_matrix(A=Mstore,Ainv=Minv,N=naos,method=1,&
    2404             :                !!!TRIM               !range1_thr=1.0E-9_dp,range2_thr=1.0E-9_dp,&
    2405             :                !!!TRIM               shift=1.0E-5_dp,&
    2406             :                !!!TRIM               range1=nmos(idomain),range2=nmos(idomain),&
    2407             :                !!!TRIM
    2408             :                !!!TRIM       ! apply the inverted matrix
    2409             :                !!!TRIM       ! RZK-warning this is only possible when the preconditioner is applied
    2410             :                !!!TRIM       tmp(:,ii)=MATMUL(Minv,subm_in(idomain)%mdata(:,ii))
    2411             :                !!!TRIM    ENDDO
    2412             :                !!!TRIM    subm_out=MATMUL(tmp,sigma)
    2413             :                !!!TRIM    deallocate(tmp)
    2414             :                !!!TRIM ELSE
    2415             : 
    2416         525 :                IF (PRESENT(bad_modes_projector_down)) THEN
    2417         240 :                   ALLOCATE (proj_array(naos, naos))
    2418             :                   CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, &
    2419             :                                             range1=nmos(idomain), range2=n_domain_mos, &
    2420             :                                             range1_thr=eps_zero_eigenvalues, &
    2421             :                                             bad_modes_projector_down=proj_array, &
    2422             :                                             s_inv_half=subm_s_inv_half(idomain)%mdata, &
    2423             :                                             s_half=subm_s_half(idomain)%mdata &
    2424          60 :                                             )
    2425             :                ELSE
    2426             :                   CALL pseudo_invert_matrix(A=subm_main(idomain)%mdata, Ainv=Minv, N=naos, method=1, &
    2427         465 :                                             range1=nmos(idomain), range2=n_domain_mos)
    2428             :                END IF
    2429             :                !!!TRIM ENDIF
    2430             : 
    2431         525 :                CALL copy_submatrices(subm_main(idomain), preconditioner(idomain), .FALSE.)
    2432         525 :                CALL copy_submatrix_data(Minv, preconditioner(idomain))
    2433         525 :                DEALLOCATE (Minv)
    2434             : 
    2435         525 :                IF (PRESENT(bad_modes_projector_down)) THEN
    2436          60 :                   CALL copy_submatrices(subm_main(idomain), bad_modes_projector_down(idomain), .FALSE.)
    2437          60 :                   CALL copy_submatrix_data(proj_array, bad_modes_projector_down(idomain))
    2438          60 :                   DEALLOCATE (proj_array)
    2439             :                END IF
    2440             : 
    2441             :             END IF ! submatrix for the domain exists
    2442             : 
    2443             :          END DO ! loop over domains
    2444             : 
    2445             :       END IF
    2446             : 
    2447         550 :       CALL release_submatrices(subm_main)
    2448        3416 :       DEALLOCATE (subm_main)
    2449             :       !DEALLOCATE(subm_s)
    2450             :       !DEALLOCATE(subm_r)
    2451             : 
    2452             :       !IF (matrix_r_required) THEN
    2453             :       !   CALL dbcsr_release(m_tmp_no_1)
    2454             :       !   CALL dbcsr_release(m_tmp_no_2)
    2455             :       !   CALL dbcsr_release(matrix_r)
    2456             :       !ENDIF
    2457             : 
    2458         550 :       CALL timestop(handle)
    2459             : 
    2460        1650 :    END SUBROUTINE construct_domain_preconditioner
    2461             : 
    2462             : ! **************************************************************************************************
    2463             : !> \brief Constructs S^(+1/2) and S^(-1/2) submatrices for each domain
    2464             : !> \param matrix_s ...
    2465             : !> \param subm_s_sqrt ...
    2466             : !> \param subm_s_sqrt_inv ...
    2467             : !> \param dpattern ...
    2468             : !> \param map ...
    2469             : !> \param node_of_domain ...
    2470             : !> \par History
    2471             : !>       2013.03 created [Rustam Z. Khaliullin]
    2472             : !> \author Rustam Z. Khaliullin
    2473             : ! **************************************************************************************************
    2474          40 :    SUBROUTINE construct_domain_s_sqrt(matrix_s, subm_s_sqrt, subm_s_sqrt_inv, &
    2475          40 :                                       dpattern, map, node_of_domain)
    2476             : 
    2477             :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_s
    2478             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2479             :          INTENT(INOUT)                                   :: subm_s_sqrt, subm_s_sqrt_inv
    2480             :       TYPE(dbcsr_type), INTENT(IN)                       :: dpattern
    2481             :       TYPE(domain_map_type), INTENT(IN)                  :: map
    2482             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: node_of_domain
    2483             : 
    2484             :       CHARACTER(len=*), PARAMETER :: routineN = 'construct_domain_s_sqrt'
    2485             : 
    2486             :       INTEGER                                            :: handle, idomain, naos, ndomains
    2487          40 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Ssqrt, Ssqrtinv
    2488             :       TYPE(domain_submatrix_type), ALLOCATABLE, &
    2489             :          DIMENSION(:)                                    :: subm_s
    2490             : 
    2491          40 :       CALL timeset(routineN, handle)
    2492             : 
    2493          40 :       ndomains = dbcsr_nblkcols_total(dpattern)
    2494          40 :       CPASSERT(SIZE(subm_s_sqrt) .EQ. ndomains)
    2495          40 :       CPASSERT(SIZE(subm_s_sqrt_inv) .EQ. ndomains)
    2496         376 :       ALLOCATE (subm_s(ndomains))
    2497          40 :       CALL init_submatrices(subm_s)
    2498             : 
    2499             :       CALL construct_submatrices(matrix_s, subm_s, &
    2500          40 :                                  dpattern, map, node_of_domain, select_row_col)
    2501             : 
    2502             :       ! loop over domains - perform inversion
    2503         296 :       DO idomain = 1, ndomains
    2504             : 
    2505             :          ! check if the submatrix exists
    2506         296 :          IF (subm_s(idomain)%domain .GT. 0) THEN
    2507             : 
    2508         128 :             naos = subm_s(idomain)%nrows
    2509             : 
    2510         512 :             ALLOCATE (Ssqrt(naos, naos))
    2511         512 :             ALLOCATE (Ssqrtinv(naos, naos))
    2512             : 
    2513             :             CALL matrix_sqrt(A=subm_s(idomain)%mdata, Asqrt=Ssqrt, Asqrtinv=Ssqrtinv, &
    2514         128 :                              N=naos)
    2515             : 
    2516         128 :             CALL copy_submatrices(subm_s(idomain), subm_s_sqrt(idomain), .FALSE.)
    2517         128 :             CALL copy_submatrix_data(Ssqrt, subm_s_sqrt(idomain))
    2518             : 
    2519         128 :             CALL copy_submatrices(subm_s(idomain), subm_s_sqrt_inv(idomain), .FALSE.)
    2520         128 :             CALL copy_submatrix_data(Ssqrtinv, subm_s_sqrt_inv(idomain))
    2521             : 
    2522         128 :             DEALLOCATE (Ssqrtinv)
    2523         128 :             DEALLOCATE (Ssqrt)
    2524             : 
    2525             :          END IF ! submatrix for the domain exists
    2526             : 
    2527             :       END DO ! loop over domains
    2528             : 
    2529          40 :       CALL release_submatrices(subm_s)
    2530         296 :       DEALLOCATE (subm_s)
    2531             : 
    2532          40 :       CALL timestop(handle)
    2533             : 
    2534          80 :    END SUBROUTINE construct_domain_s_sqrt
    2535             : 
    2536             : ! **************************************************************************************************
    2537             : !> \brief Constructs S_inv block for each domain
    2538             : !> \param matrix_s ...
    2539             : !> \param subm_s_inv ...
    2540             : !> \param dpattern ...
    2541             : !> \param map ...
    2542             : !> \param node_of_domain ...
    2543             : !> \par History
    2544             : !>       2013.02 created [Rustam Z. Khaliullin]
    2545             : !> \author Rustam Z. Khaliullin
    2546             : ! **************************************************************************************************
    2547          54 :    SUBROUTINE construct_domain_s_inv(matrix_s, subm_s_inv, dpattern, map, &
    2548          54 :                                      node_of_domain)
    2549             :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_s
    2550             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2551             :          INTENT(INOUT)                                   :: subm_s_inv
    2552             :       TYPE(dbcsr_type), INTENT(IN)                       :: dpattern
    2553             :       TYPE(domain_map_type), INTENT(IN)                  :: map
    2554             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: node_of_domain
    2555             : 
    2556             :       CHARACTER(len=*), PARAMETER :: routineN = 'construct_domain_s_inv'
    2557             : 
    2558             :       INTEGER                                            :: handle, idomain, naos, ndomains
    2559          54 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Sinv
    2560             :       TYPE(domain_submatrix_type), ALLOCATABLE, &
    2561             :          DIMENSION(:)                                    :: subm_s
    2562             : 
    2563          54 :       CALL timeset(routineN, handle)
    2564             : 
    2565          54 :       ndomains = dbcsr_nblkcols_total(dpattern)
    2566             : 
    2567          54 :       CPASSERT(SIZE(subm_s_inv) .EQ. ndomains)
    2568         520 :       ALLOCATE (subm_s(ndomains))
    2569          54 :       CALL init_submatrices(subm_s)
    2570             : 
    2571             :       CALL construct_submatrices(matrix_s, subm_s, &
    2572          54 :                                  dpattern, map, node_of_domain, select_row_col)
    2573             : 
    2574             :       ! loop over domains - perform inversion
    2575         412 :       DO idomain = 1, ndomains
    2576             : 
    2577             :          ! check if the submatrix exists
    2578         412 :          IF (subm_s(idomain)%domain .GT. 0) THEN
    2579             : 
    2580         179 :             naos = subm_s(idomain)%nrows
    2581             : 
    2582         716 :             ALLOCATE (Sinv(naos, naos))
    2583             : 
    2584             :             CALL pseudo_invert_matrix(A=subm_s(idomain)%mdata, Ainv=Sinv, N=naos, &
    2585         179 :                                       method=0)
    2586             : 
    2587         179 :             CALL copy_submatrices(subm_s(idomain), subm_s_inv(idomain), .FALSE.)
    2588         179 :             CALL copy_submatrix_data(Sinv, subm_s_inv(idomain))
    2589             : 
    2590         179 :             DEALLOCATE (Sinv)
    2591             : 
    2592             :          END IF ! submatrix for the domain exists
    2593             : 
    2594             :       END DO ! loop over domains
    2595             : 
    2596          54 :       CALL release_submatrices(subm_s)
    2597         412 :       DEALLOCATE (subm_s)
    2598             : 
    2599          54 :       CALL timestop(handle)
    2600             : 
    2601         108 :    END SUBROUTINE construct_domain_s_inv
    2602             : 
    2603             : ! **************************************************************************************************
    2604             : !> \brief Constructs subblocks of the covariant-covariant projectors (i.e. DM without spin factor)
    2605             : !> \param matrix_t ...
    2606             : !> \param matrix_sigma_inv ...
    2607             : !> \param matrix_s ...
    2608             : !> \param subm_r_down ...
    2609             : !> \param dpattern ...
    2610             : !> \param map ...
    2611             : !> \param node_of_domain ...
    2612             : !> \param filter_eps ...
    2613             : !> \par History
    2614             : !>       2013.02 created [Rustam Z. Khaliullin]
    2615             : !> \author Rustam Z. Khaliullin
    2616             : ! **************************************************************************************************
    2617          48 :    SUBROUTINE construct_domain_r_down(matrix_t, matrix_sigma_inv, matrix_s, &
    2618          24 :                                       subm_r_down, dpattern, map, node_of_domain, filter_eps)
    2619             : 
    2620             :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_t, matrix_sigma_inv, matrix_s
    2621             :       TYPE(domain_submatrix_type), DIMENSION(:), &
    2622             :          INTENT(INOUT)                                   :: subm_r_down
    2623             :       TYPE(dbcsr_type), INTENT(IN)                       :: dpattern
    2624             :       TYPE(domain_map_type), INTENT(IN)                  :: map
    2625             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: node_of_domain
    2626             :       REAL(KIND=dp)                                      :: filter_eps
    2627             : 
    2628             :       CHARACTER(len=*), PARAMETER :: routineN = 'construct_domain_r_down'
    2629             : 
    2630             :       INTEGER                                            :: handle, ndomains
    2631             :       TYPE(dbcsr_type)                                   :: m_tmp_no_1, m_tmp_no_2, matrix_r
    2632             : 
    2633          24 :       CALL timeset(routineN, handle)
    2634             : 
    2635             :       ! compute the density matrix in the COVARIANT representation
    2636             :       CALL dbcsr_create(matrix_r, &
    2637             :                         template=matrix_s, &
    2638          24 :                         matrix_type=dbcsr_type_symmetric)
    2639             :       CALL dbcsr_create(m_tmp_no_1, &
    2640             :                         template=matrix_t, &
    2641          24 :                         matrix_type=dbcsr_type_no_symmetry)
    2642             :       CALL dbcsr_create(m_tmp_no_2, &
    2643             :                         template=matrix_t, &
    2644          24 :                         matrix_type=dbcsr_type_no_symmetry)
    2645             : 
    2646             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s, matrix_t, &
    2647          24 :                           0.0_dp, m_tmp_no_1, filter_eps=filter_eps)
    2648             :       CALL dbcsr_multiply("N", "N", 1.0_dp, m_tmp_no_1, matrix_sigma_inv, &
    2649          24 :                           0.0_dp, m_tmp_no_2, filter_eps=filter_eps)
    2650             :       CALL dbcsr_multiply("N", "T", 1.0_dp, m_tmp_no_2, m_tmp_no_1, &
    2651          24 :                           0.0_dp, matrix_r, filter_eps=filter_eps)
    2652             : 
    2653          24 :       CALL dbcsr_release(m_tmp_no_1)
    2654          24 :       CALL dbcsr_release(m_tmp_no_2)
    2655             : 
    2656          24 :       ndomains = dbcsr_nblkcols_total(dpattern)
    2657          24 :       CPASSERT(SIZE(subm_r_down) .EQ. ndomains)
    2658             : 
    2659             :       CALL construct_submatrices(matrix_r, subm_r_down, &
    2660          24 :                                  dpattern, map, node_of_domain, select_row_col)
    2661             : 
    2662          24 :       CALL dbcsr_release(matrix_r)
    2663             : 
    2664          24 :       CALL timestop(handle)
    2665             : 
    2666          24 :    END SUBROUTINE construct_domain_r_down
    2667             : 
    2668             : ! **************************************************************************************************
    2669             : !> \brief Finds the square root of a matrix and its inverse
    2670             : !> \param A ...
    2671             : !> \param Asqrt ...
    2672             : !> \param Asqrtinv ...
    2673             : !> \param N ...
    2674             : !> \par History
    2675             : !>       2013.03 created [Rustam Z. Khaliullin]
    2676             : !> \author Rustam Z. Khaliullin
    2677             : ! **************************************************************************************************
    2678         128 :    SUBROUTINE matrix_sqrt(A, Asqrt, Asqrtinv, N)
    2679             : 
    2680             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: A
    2681             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: Asqrt, Asqrtinv
    2682             :       INTEGER, INTENT(IN)                                :: N
    2683             : 
    2684             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'matrix_sqrt'
    2685             : 
    2686             :       INTEGER                                            :: handle, INFO, jj, LWORK, unit_nr
    2687         128 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, WORK
    2688         128 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: test, testN
    2689             :       TYPE(cp_logger_type), POINTER                      :: logger
    2690             : 
    2691         128 :       CALL timeset(routineN, handle)
    2692             : 
    2693      585400 :       Asqrtinv = A
    2694         128 :       INFO = 0
    2695             : 
    2696             :       ! get a useful unit_nr
    2697         128 :       logger => cp_get_default_logger()
    2698         128 :       IF (logger%para_env%is_source()) THEN
    2699          65 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    2700             :       ELSE
    2701             :          unit_nr = -1
    2702             :       END IF
    2703             : 
    2704             :       !CALL DPOTRF('L', N, Ainv, N, INFO )
    2705             :       !IF( INFO.NE.0 ) THEN
    2706             :       !   CPErrorMessage(cp_failure_level,routineP,"DPOTRF failed")
    2707             :       !   CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
    2708             :       !END IF
    2709             :       !CALL DPOTRI('L', N, Ainv, N, INFO )
    2710             :       !IF( INFO.NE.0 ) THEN
    2711             :       !   CPErrorMessage(cp_failure_level,routineP,"DPOTRI failed")
    2712             :       !   CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
    2713             :       !END IF
    2714             :       !! complete the matrix
    2715             :       !DO ii=1,N
    2716             :       !   DO jj=ii+1,N
    2717             :       !      Ainv(ii,jj)=Ainv(jj,ii)
    2718             :       !   ENDDO
    2719             :       !   !WRITE(*,'(100F13.9)') Ainv(ii,:)
    2720             :       !ENDDO
    2721             : 
    2722             :       ! diagonalize first
    2723         384 :       ALLOCATE (eigenvalues(N))
    2724             :       ! Query the optimal workspace for dsyev
    2725         128 :       LWORK = -1
    2726         128 :       ALLOCATE (WORK(MAX(1, LWORK)))
    2727         128 :       CALL DSYEV('V', 'L', N, Asqrtinv, N, eigenvalues, WORK, LWORK, INFO)
    2728         128 :       LWORK = INT(WORK(1))
    2729         128 :       DEALLOCATE (WORK)
    2730             :       ! Allocate the workspace and solve the eigenproblem
    2731         384 :       ALLOCATE (WORK(MAX(1, LWORK)))
    2732         128 :       CALL DSYEV('V', 'L', N, Asqrtinv, N, eigenvalues, WORK, LWORK, INFO)
    2733         128 :       IF (INFO .NE. 0) THEN
    2734           0 :          IF (unit_nr > 0) WRITE (unit_nr, *) 'DSYEV ERROR MESSAGE: ', INFO
    2735           0 :          CPABORT("DSYEV failed")
    2736             :       END IF
    2737         128 :       DEALLOCATE (WORK)
    2738             : 
    2739             :       ! take functions of eigenvalues and use eigenvectors to compute the matrix function
    2740             :       ! first sqrt
    2741         512 :       ALLOCATE (test(N, N))
    2742        6242 :       DO jj = 1, N
    2743      585400 :          test(jj, :) = Asqrtinv(:, jj)*SQRT(eigenvalues(jj))
    2744             :       END DO
    2745         512 :       ALLOCATE (testN(N, N))
    2746      899950 :       testN(:, :) = MATMUL(Asqrtinv, test)
    2747      585400 :       Asqrt = testN
    2748             :       ! now, sqrt_inv
    2749        6242 :       DO jj = 1, N
    2750      585400 :          test(jj, :) = Asqrtinv(:, jj)/SQRT(eigenvalues(jj))
    2751             :       END DO
    2752      899950 :       testN(:, :) = MATMUL(Asqrtinv, test)
    2753      585400 :       Asqrtinv = testN
    2754         128 :       DEALLOCATE (test, testN)
    2755             : 
    2756         128 :       DEALLOCATE (eigenvalues)
    2757             : 
    2758             : !!!    ! compute the error
    2759             : !!!    allocate(test(N,N))
    2760             : !!!    test=MATMUL(Ainv,A)
    2761             : !!!    DO ii=1,N
    2762             : !!!       test(ii,ii)=test(ii,ii)-1.0_dp
    2763             : !!!    ENDDO
    2764             : !!!    test_error=0.0_dp
    2765             : !!!    DO ii=1,N
    2766             : !!!       DO jj=1,N
    2767             : !!!          test_error=test_error+test(jj,ii)*test(jj,ii)
    2768             : !!!       ENDDO
    2769             : !!!    ENDDO
    2770             : !!!    WRITE(*,*) "Inversion error: ", SQRT(test_error)
    2771             : !!!    deallocate(test)
    2772             : 
    2773         128 :       CALL timestop(handle)
    2774             : 
    2775         128 :    END SUBROUTINE matrix_sqrt
    2776             : 
    2777             : ! **************************************************************************************************
    2778             : !> \brief Inverts a matrix using a requested method
    2779             : !>         0. Cholesky factorization
    2780             : !>         1. Diagonalization
    2781             : !> \param A ...
    2782             : !> \param Ainv ...
    2783             : !> \param N ...
    2784             : !> \param method ...
    2785             : !> \param range1 ...
    2786             : !> \param range2 ...
    2787             : !> \param range1_thr ...
    2788             : !> \param shift ...
    2789             : !> \param bad_modes_projector_down ...
    2790             : !> \param s_inv_half ...
    2791             : !> \param s_half ...
    2792             : !> \par History
    2793             : !>       2012.04 created [Rustam Z. Khaliullin]
    2794             : !> \author Rustam Z. Khaliullin
    2795             : ! **************************************************************************************************
    2796        1023 :    SUBROUTINE pseudo_invert_matrix(A, Ainv, N, method, range1, range2, range1_thr, &
    2797        2046 :                                    shift, bad_modes_projector_down, s_inv_half, s_half)
    2798             : 
    2799             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: A
    2800             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: Ainv
    2801             :       INTEGER, INTENT(IN)                                :: N, method
    2802             :       INTEGER, INTENT(IN), OPTIONAL                      :: range1, range2
    2803             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: range1_thr, shift
    2804             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
    2805             :          OPTIONAL                                        :: bad_modes_projector_down
    2806             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
    2807             :          OPTIONAL                                        :: s_inv_half, s_half
    2808             : 
    2809             :       CHARACTER(len=*), PARAMETER :: routineN = 'pseudo_invert_matrix'
    2810             : 
    2811             :       INTEGER                                            :: handle, ii, INFO, jj, LWORK, range1_eiv, &
    2812             :                                                             range2_eiv, range3_eiv, unit_nr
    2813             :       LOGICAL                                            :: use_both, use_ranges_only, use_thr_only
    2814             :       REAL(KIND=dp)                                      :: my_shift
    2815        1023 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, WORK
    2816        1023 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: temp1, temp2, temp3, temp4
    2817             :       TYPE(cp_logger_type), POINTER                      :: logger
    2818             : 
    2819        1023 :       CALL timeset(routineN, handle)
    2820             : 
    2821             :       ! get a useful unit_nr
    2822        1023 :       logger => cp_get_default_logger()
    2823        1023 :       IF (logger%para_env%is_source()) THEN
    2824         514 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    2825             :       ELSE
    2826             :          unit_nr = -1
    2827             :       END IF
    2828             : 
    2829        1023 :       IF (method .EQ. 1) THEN
    2830             : 
    2831         844 :          IF ((PRESENT(range1) .AND. (.NOT. PRESENT(range2))) .OR. (PRESENT(range2) .AND. (.NOT. PRESENT(range1)))) THEN
    2832           0 :             CPABORT("range1 and range2 must be provided together")
    2833             :          END IF
    2834             : 
    2835         844 :          IF (PRESENT(range1) .AND. PRESENT(range1_thr)) THEN
    2836             :             use_both = .TRUE.
    2837             :             use_thr_only = .FALSE.
    2838             :             use_ranges_only = .FALSE.
    2839             :          ELSE
    2840         784 :             use_both = .FALSE.
    2841             : 
    2842         784 :             IF (PRESENT(range1)) THEN
    2843             :                use_ranges_only = .TRUE.
    2844             :             ELSE
    2845           0 :                use_ranges_only = .FALSE.
    2846             :             END IF
    2847             : 
    2848         784 :             IF (PRESENT(range1_thr)) THEN
    2849             :                use_thr_only = .TRUE.
    2850             :             ELSE
    2851         784 :                use_thr_only = .FALSE.
    2852             :             END IF
    2853             : 
    2854             :          END IF
    2855             : 
    2856         844 :          IF ((PRESENT(s_half) .AND. (.NOT. PRESENT(s_inv_half))) .OR. (PRESENT(s_inv_half) .AND. (.NOT. PRESENT(s_half)))) THEN
    2857           0 :             CPABORT("Domain overlap matrix missing")
    2858             :          END IF
    2859             :       END IF
    2860             : 
    2861        1023 :       my_shift = 0.0_dp
    2862        1023 :       IF (PRESENT(shift)) THEN
    2863         319 :          my_shift = shift
    2864             :       END IF
    2865             : 
    2866     2592719 :       Ainv = A
    2867        1023 :       INFO = 0
    2868             : 
    2869         179 :       SELECT CASE (method)
    2870             :       CASE (0) ! Inversion via cholesky factorization
    2871             : 
    2872         179 :          CALL DPOTRF('L', N, Ainv, N, INFO)
    2873         179 :          IF (INFO .NE. 0) THEN
    2874           0 :             CPABORT("DPOTRF failed")
    2875             :          END IF
    2876         179 :          CALL DPOTRI('L', N, Ainv, N, INFO)
    2877         179 :          IF (INFO .NE. 0) THEN
    2878           0 :             CPABORT("DPOTRI failed")
    2879             :          END IF
    2880             :          ! complete the matrix
    2881        6834 :          DO ii = 1, N
    2882      288702 :             DO jj = ii + 1, N
    2883      288523 :                Ainv(ii, jj) = Ainv(jj, ii)
    2884             :             END DO
    2885             :             !WRITE(*,'(100F13.9)') Ainv(ii,:)
    2886             :          END DO
    2887             : 
    2888             :       CASE (1)
    2889             : 
    2890             :          ! diagonalize first
    2891        2532 :          ALLOCATE (eigenvalues(N))
    2892        3376 :          ALLOCATE (temp1(N, N))
    2893        3376 :          ALLOCATE (temp4(N, N))
    2894         844 :          IF (PRESENT(s_inv_half)) THEN
    2895          60 :             CALL DSYMM('L', 'U', N, N, 1.0_dp, s_inv_half, N, A, N, 0.0_dp, temp1, N)
    2896          60 :             CALL DSYMM('R', 'U', N, N, 1.0_dp, s_inv_half, N, temp1, N, 0.0_dp, Ainv, N)
    2897             :          END IF
    2898             :          ! Query the optimal workspace for dsyev
    2899         844 :          LWORK = -1
    2900         844 :          ALLOCATE (WORK(MAX(1, LWORK)))
    2901         844 :          CALL DSYEV('V', 'L', N, Ainv, N, eigenvalues, WORK, LWORK, INFO)
    2902             : 
    2903         844 :          LWORK = INT(WORK(1))
    2904         844 :          DEALLOCATE (WORK)
    2905             :          ! Allocate the workspace and solve the eigenproblem
    2906        2532 :          ALLOCATE (WORK(MAX(1, LWORK)))
    2907         844 :          CALL DSYEV('V', 'L', N, Ainv, N, eigenvalues, WORK, LWORK, INFO)
    2908             : 
    2909         844 :          IF (INFO .NE. 0) THEN
    2910           0 :             IF (unit_nr > 0) WRITE (unit_nr, *) 'EIGENSYSTEM ERROR MESSAGE: ', INFO
    2911           0 :             CPABORT("Eigenproblem routine failed")
    2912             :          END IF
    2913         844 :          DEALLOCATE (WORK)
    2914             : 
    2915             :          !WRITE(*,*) "EIGENVALS: "
    2916             :          !WRITE(*,'(4F13.9)') eigenvalues(:)
    2917             : 
    2918             :          ! invert eigenvalues and use eigenvectors to compute pseudo Ainv
    2919             :          ! project out near-zero eigenvalue modes
    2920        3376 :          ALLOCATE (temp2(N, N))
    2921        1024 :          IF (PRESENT(bad_modes_projector_down)) ALLOCATE (temp3(N, N))
    2922     2015494 :          temp2(1:N, 1:N) = Ainv(1:N, 1:N)
    2923             : 
    2924         844 :          range1_eiv = 0
    2925         844 :          range2_eiv = 0
    2926         844 :          range3_eiv = 0
    2927             : 
    2928         844 :          IF (use_both) THEN
    2929        1116 :             DO jj = 1, N
    2930        1116 :                IF ((jj .LE. range2) .AND. (eigenvalues(jj) .LT. range1_thr)) THEN
    2931        9274 :                   temp1(jj, :) = temp2(:, jj)*0.0_dp
    2932        9274 :                   IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp
    2933         454 :                   range1_eiv = range1_eiv + 1
    2934             :                ELSE
    2935       17528 :                   temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
    2936       17528 :                   IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp
    2937         602 :                   range2_eiv = range2_eiv + 1
    2938             :                END IF
    2939             :             END DO
    2940             :          ELSE
    2941         784 :             IF (use_ranges_only) THEN
    2942       34678 :                DO jj = 1, N
    2943       34678 :                   IF (jj .LE. range1) THEN
    2944      152199 :                      temp1(jj, :) = temp2(:, jj)*0.0_dp
    2945        3173 :                      IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp
    2946        3173 :                      range1_eiv = range1_eiv + 1
    2947       30721 :                   ELSE IF (jj .LE. range2) THEN
    2948      256093 :                      temp1(jj, :) = temp2(:, jj)*1.0_dp
    2949        4129 :                      IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp
    2950        4129 :                      range2_eiv = range2_eiv + 1
    2951             :                   ELSE
    2952     1579556 :                      temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
    2953       26592 :                      IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp
    2954       26592 :                      range3_eiv = range3_eiv + 1
    2955             :                   END IF
    2956             :                END DO
    2957           0 :             ELSE IF (use_thr_only) THEN
    2958           0 :                DO jj = 1, N
    2959           0 :                   IF (eigenvalues(jj) .LT. range1_thr) THEN
    2960           0 :                      temp1(jj, :) = temp2(:, jj)*0.0_dp
    2961           0 :                      IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*1.0_dp
    2962           0 :                      range1_eiv = range1_eiv + 1
    2963             :                   ELSE
    2964           0 :                      temp1(jj, :) = temp2(:, jj)/(eigenvalues(jj) + my_shift)
    2965           0 :                      IF (PRESENT(bad_modes_projector_down)) temp3(jj, :) = Ainv(:, jj)*0.0_dp
    2966           0 :                      range2_eiv = range2_eiv + 1
    2967             :                   END IF
    2968             :                END DO
    2969             :             ELSE ! no ranges, no thresholds
    2970           0 :                CPABORT("Invert using Cholesky. It would be faster.")
    2971             :             END IF
    2972             :          END IF
    2973             :          !WRITE(*,*) ' EIV RANGES: ', range1_eiv, range2_eiv, range3_eiv
    2974         844 :          IF (PRESENT(bad_modes_projector_down)) THEN
    2975          60 :             IF (PRESENT(s_half)) THEN
    2976          60 :                CALL DSYMM('L', 'U', N, N, 1.0_dp, s_half, N, temp2, N, 0.0_dp, Ainv, N)
    2977          60 :                CALL DSYMM('R', 'U', N, N, 1.0_dp, s_half, N, temp3, N, 0.0_dp, temp4, N)
    2978          60 :                CALL DGEMM('N', 'N', N, N, N, 1.0_dp, Ainv, N, temp4, N, 0.0_dp, bad_modes_projector_down, N)
    2979             :             ELSE
    2980           0 :                CALL DGEMM('N', 'N', N, N, N, 1.0_dp, temp2, N, temp3, N, 0.0_dp, bad_modes_projector_down, N)
    2981             :             END IF
    2982             :          END IF
    2983             : 
    2984         844 :          IF (PRESENT(s_inv_half)) THEN
    2985          60 :             CALL DSYMM('L', 'U', N, N, 1.0_dp, s_inv_half, N, temp2, N, 0.0_dp, temp4, N)
    2986          60 :             CALL DSYMM('R', 'U', N, N, 1.0_dp, s_inv_half, N, temp1, N, 0.0_dp, temp2, N)
    2987          60 :             CALL DGEMM('N', 'N', N, N, N, 1.0_dp, temp4, N, temp2, N, 0.0_dp, Ainv, N)
    2988             :          ELSE
    2989         784 :             CALL DGEMM('N', 'N', N, N, N, 1.0_dp, temp2, N, temp1, N, 0.0_dp, Ainv, N)
    2990             :          END IF
    2991         844 :          DEALLOCATE (temp1, temp2, temp4)
    2992         844 :          IF (PRESENT(bad_modes_projector_down)) DEALLOCATE (temp3)
    2993         844 :          DEALLOCATE (eigenvalues)
    2994             : 
    2995             :       CASE DEFAULT
    2996             : 
    2997        1023 :          CPABORT("Illegal method selected for matrix inversion")
    2998             : 
    2999             :       END SELECT
    3000             : 
    3001             :       !! compute the inversion error
    3002             :       !allocate(temp1(N,N))
    3003             :       !temp1=MATMUL(Ainv,A)
    3004             :       !DO ii=1,N
    3005             :       !   temp1(ii,ii)=temp1(ii,ii)-1.0_dp
    3006             :       !ENDDO
    3007             :       !temp1_error=0.0_dp
    3008             :       !DO ii=1,N
    3009             :       !   DO jj=1,N
    3010             :       !      temp1_error=temp1_error+temp1(jj,ii)*temp1(jj,ii)
    3011             :       !   ENDDO
    3012             :       !ENDDO
    3013             :       !WRITE(*,*) "Inversion error: ", SQRT(temp1_error)
    3014             :       !deallocate(temp1)
    3015             : 
    3016        1023 :       CALL timestop(handle)
    3017             : 
    3018        1023 :    END SUBROUTINE pseudo_invert_matrix
    3019             : 
    3020             : ! **************************************************************************************************
    3021             : !> \brief Find matrix power using diagonalization
    3022             : !> \param A ...
    3023             : !> \param Apow ...
    3024             : !> \param power ...
    3025             : !> \param N ...
    3026             : !> \param range1 ...
    3027             : !> \param range1_thr ...
    3028             : !> \param shift ...
    3029             : !> \par History
    3030             : !>       2012.04 created [Rustam Z. Khaliullin]
    3031             : !> \author Rustam Z. Khaliullin
    3032             : ! **************************************************************************************************
    3033             : 
    3034           0 :    SUBROUTINE pseudo_matrix_power(A, Apow, power, N, range1, range1_thr, shift)
    3035             : 
    3036             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: A
    3037             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: Apow
    3038             :       REAL(KIND=dp), INTENT(IN)                          :: power
    3039             :       INTEGER, INTENT(IN)                                :: N
    3040             :       INTEGER, INTENT(IN), OPTIONAL                      :: range1
    3041             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: range1_thr, shift
    3042             : 
    3043             :       CHARACTER(len=*), PARAMETER :: routineN = 'pseudo_matrix_power'
    3044             : 
    3045             :       INTEGER                                            :: handle, INFO, jj, LWORK, range1_eiv, &
    3046             :                                                             range2_eiv, unit_nr
    3047             :       LOGICAL                                            :: use_both, use_ranges, use_thr
    3048             :       REAL(KIND=dp)                                      :: my_shift
    3049           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, WORK
    3050           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: temp1, temp2
    3051             :       TYPE(cp_logger_type), POINTER                      :: logger
    3052             : 
    3053           0 :       CALL timeset(routineN, handle)
    3054             : 
    3055             :       ! get a useful unit_nr
    3056           0 :       logger => cp_get_default_logger()
    3057           0 :       IF (logger%para_env%is_source()) THEN
    3058           0 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    3059             :       ELSE
    3060             :          unit_nr = -1
    3061             :       END IF
    3062             : 
    3063           0 :       IF (PRESENT(range1) .AND. PRESENT(range1_thr)) THEN
    3064             :          use_both = .TRUE.
    3065             :       ELSE
    3066           0 :          use_both = .FALSE.
    3067           0 :          IF (PRESENT(range1)) THEN
    3068             :             use_ranges = .TRUE.
    3069             :          ELSE
    3070           0 :             use_ranges = .FALSE.
    3071           0 :             IF (PRESENT(range1_thr)) THEN
    3072             :                use_thr = .TRUE.
    3073             :             ELSE
    3074           0 :                use_thr = .FALSE.
    3075             :             END IF
    3076             :          END IF
    3077             :       END IF
    3078             : 
    3079           0 :       my_shift = 0.0_dp
    3080           0 :       IF (PRESENT(shift)) THEN
    3081           0 :          my_shift = shift
    3082             :       END IF
    3083             : 
    3084           0 :       Apow = A
    3085           0 :       INFO = 0
    3086             : 
    3087             :       ! diagonalize first
    3088           0 :       ALLOCATE (eigenvalues(N))
    3089           0 :       ALLOCATE (temp1(N, N))
    3090             : 
    3091             :       ! Query the optimal workspace for dsyev
    3092           0 :       LWORK = -1
    3093           0 :       ALLOCATE (WORK(MAX(1, LWORK)))
    3094           0 :       CALL DSYEV('V', 'L', N, Apow, N, eigenvalues, WORK, LWORK, INFO)
    3095             : 
    3096           0 :       LWORK = INT(WORK(1))
    3097           0 :       DEALLOCATE (WORK)
    3098             :       ! Allocate the workspace and solve the eigenproblem
    3099           0 :       ALLOCATE (WORK(MAX(1, LWORK)))
    3100             : 
    3101           0 :       CALL DSYEV('V', 'L', N, Apow, N, eigenvalues, WORK, LWORK, INFO)
    3102             : 
    3103           0 :       IF (INFO .NE. 0) THEN
    3104           0 :          IF (unit_nr > 0) WRITE (unit_nr, *) 'EIGENSYSTEM ERROR MESSAGE: ', INFO
    3105           0 :          CPABORT("Eigenproblem routine failed")
    3106             :       END IF
    3107           0 :       DEALLOCATE (WORK)
    3108             : 
    3109             :       !WRITE(*,*) "EIGENVALS: "
    3110             :       !WRITE(*,'(4F13.9)') eigenvalues(:)
    3111             : 
    3112             :       ! invert eigenvalues and use eigenvectors to compute pseudo Ainv
    3113             :       ! project out near-zero eigenvalue modes
    3114           0 :       ALLOCATE (temp2(N, N))
    3115             : 
    3116           0 :       temp2(1:N, 1:N) = Apow(1:N, 1:N)
    3117             : 
    3118           0 :       range1_eiv = 0
    3119           0 :       range2_eiv = 0
    3120             : 
    3121           0 :       IF (use_both) THEN
    3122           0 :          DO jj = 1, N
    3123           0 :             IF ((jj .LE. range1) .AND. (eigenvalues(jj) .LT. range1_thr)) THEN
    3124           0 :                temp1(jj, :) = temp2(:, jj)*0.0_dp
    3125           0 :                range1_eiv = range1_eiv + 1
    3126             :             ELSE
    3127           0 :                temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
    3128             :             END IF
    3129             :          END DO
    3130             :       ELSE
    3131           0 :          IF (use_ranges) THEN
    3132           0 :             DO jj = 1, N
    3133           0 :                IF (jj .LE. range1) THEN
    3134           0 :                   temp1(jj, :) = temp2(:, jj)*0.0_dp
    3135           0 :                   range1_eiv = range1_eiv + 1
    3136             :                ELSE
    3137           0 :                   temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
    3138             :                END IF
    3139             :             END DO
    3140             :          ELSE
    3141           0 :             IF (use_thr) THEN
    3142           0 :                DO jj = 1, N
    3143           0 :                   IF (eigenvalues(jj) .LT. range1_thr) THEN
    3144           0 :                      temp1(jj, :) = temp2(:, jj)*0.0_dp
    3145             : 
    3146           0 :                      range1_eiv = range1_eiv + 1
    3147             :                   ELSE
    3148           0 :                      temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
    3149             :                   END IF
    3150             :                END DO
    3151             :             ELSE
    3152           0 :                DO jj = 1, N
    3153           0 :                   temp1(jj, :) = temp2(:, jj)*((eigenvalues(jj) + my_shift)**power)
    3154             :                END DO
    3155             :             END IF
    3156             :          END IF
    3157             :       END IF
    3158             :       !WRITE(*,*) ' EIV RANGES: ', range1_eiv, range2_eiv, range3_eiv
    3159           0 :       Apow = MATMUL(temp2, temp1)
    3160           0 :       DEALLOCATE (temp1, temp2)
    3161           0 :       DEALLOCATE (eigenvalues)
    3162             : 
    3163           0 :       CALL timestop(handle)
    3164             : 
    3165           0 :    END SUBROUTINE pseudo_matrix_power
    3166             : 
    3167             : ! **************************************************************************************************
    3168             : !> \brief Load balancing of the submatrix computations
    3169             : !> \param almo_scf_env ...
    3170             : !> \par History
    3171             : !>       2013.02 created [Rustam Z. Khaliullin]
    3172             : !> \author Rustam Z. Khaliullin
    3173             : ! **************************************************************************************************
    3174         116 :    SUBROUTINE distribute_domains(almo_scf_env)
    3175             : 
    3176             :       TYPE(almo_scf_env_type), INTENT(INOUT)             :: almo_scf_env
    3177             : 
    3178             :       CHARACTER(len=*), PARAMETER :: routineN = 'distribute_domains'
    3179             : 
    3180             :       INTEGER                                            :: handle, idomain, least_loaded, nao, &
    3181             :                                                             ncpus, ndomains
    3182         116 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: index0
    3183         116 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cpu_load, domain_load
    3184             :       TYPE(dbcsr_distribution_type)                      :: dist
    3185             : 
    3186         116 :       CALL timeset(routineN, handle)
    3187             : 
    3188         116 :       ndomains = almo_scf_env%ndomains
    3189         116 :       CALL dbcsr_get_info(almo_scf_env%matrix_s(1), distribution=dist)
    3190         116 :       CALL dbcsr_distribution_get(dist, numnodes=ncpus)
    3191             : 
    3192         348 :       ALLOCATE (domain_load(ndomains))
    3193         926 :       DO idomain = 1, ndomains
    3194         810 :          nao = almo_scf_env%nbasis_of_domain(idomain)
    3195         926 :          domain_load(idomain) = (nao*nao*nao)*1.0_dp
    3196             :       END DO
    3197             : 
    3198         348 :       ALLOCATE (index0(ndomains))
    3199             : 
    3200         116 :       CALL sort(domain_load, ndomains, index0)
    3201             : 
    3202         348 :       ALLOCATE (cpu_load(ncpus))
    3203         348 :       cpu_load(:) = 0.0_dp
    3204             : 
    3205         926 :       DO idomain = 1, ndomains
    3206        3240 :          least_loaded = MINLOC(cpu_load, 1)
    3207         810 :          cpu_load(least_loaded) = cpu_load(least_loaded) + domain_load(idomain)
    3208         926 :          almo_scf_env%cpu_of_domain(index0(idomain)) = least_loaded - 1
    3209             :       END DO
    3210             : 
    3211         116 :       DEALLOCATE (cpu_load)
    3212         116 :       DEALLOCATE (index0)
    3213         116 :       DEALLOCATE (domain_load)
    3214             : 
    3215         116 :       CALL timestop(handle)
    3216             : 
    3217         232 :    END SUBROUTINE distribute_domains
    3218             : 
    3219             : ! **************************************************************************************************
    3220             : !> \brief Tests construction and release of domain submatrices
    3221             : !> \param matrix_no ...
    3222             : !> \param dpattern ...
    3223             : !> \param map ...
    3224             : !> \param node_of_domain ...
    3225             : !> \par History
    3226             : !>       2013.01 created [Rustam Z. Khaliullin]
    3227             : !> \author Rustam Z. Khaliullin
    3228             : ! **************************************************************************************************
    3229           0 :    SUBROUTINE construct_test(matrix_no, dpattern, map, node_of_domain)
    3230             : 
    3231             :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_no, dpattern
    3232             :       TYPE(domain_map_type), INTENT(IN)                  :: map
    3233             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: node_of_domain
    3234             : 
    3235             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'construct_test'
    3236             : 
    3237             :       INTEGER                                            :: GroupID, handle, ndomains
    3238             :       TYPE(dbcsr_type)                                   :: copy1
    3239             :       TYPE(domain_submatrix_type), ALLOCATABLE, &
    3240           0 :          DIMENSION(:)                                    :: subm_nn, subm_no
    3241             :       TYPE(mp_comm_type)                                 :: group
    3242             : 
    3243           0 :       CALL timeset(routineN, handle)
    3244             : 
    3245           0 :       ndomains = dbcsr_nblkcols_total(dpattern)
    3246           0 :       CALL dbcsr_get_info(dpattern, group=GroupID)
    3247           0 :       CALL group%set_handle(groupid)
    3248             : 
    3249           0 :       ALLOCATE (subm_no(ndomains), subm_nn(ndomains))
    3250           0 :       CALL init_submatrices(subm_no)
    3251           0 :       CALL init_submatrices(subm_nn)
    3252             : 
    3253             :       !CALL dbcsr_print(matrix_nn)
    3254             :       !CALL construct_submatrices(matrix_nn,subm_nn,dpattern,map,select_row_col)
    3255             :       !CALL print_submatrices(subm_nn,Group)
    3256             : 
    3257             :       !CALL dbcsr_print(matrix_no)
    3258           0 :       CALL construct_submatrices(matrix_no, subm_no, dpattern, map, node_of_domain, select_row)
    3259           0 :       CALL print_submatrices(subm_no, group)
    3260             : 
    3261           0 :       CALL dbcsr_create(copy1, template=matrix_no)
    3262           0 :       CALL dbcsr_copy(copy1, matrix_no)
    3263           0 :       CALL dbcsr_print(copy1)
    3264           0 :       CALL construct_dbcsr_from_submatrices(copy1, subm_no, dpattern)
    3265           0 :       CALL dbcsr_print(copy1)
    3266           0 :       CALL dbcsr_release(copy1)
    3267             : 
    3268           0 :       CALL release_submatrices(subm_no)
    3269           0 :       CALL release_submatrices(subm_nn)
    3270           0 :       DEALLOCATE (subm_no, subm_nn)
    3271             : 
    3272           0 :       CALL timestop(handle)
    3273             : 
    3274           0 :    END SUBROUTINE construct_test
    3275             : 
    3276             : ! **************************************************************************************************
    3277             : !> \brief create the initial guess for XALMOs
    3278             : !> \param m_guess ...
    3279             : !> \param m_t_in ...
    3280             : !> \param m_t0 ...
    3281             : !> \param m_quench_t ...
    3282             : !> \param m_overlap ...
    3283             : !> \param m_sigma_tmpl ...
    3284             : !> \param nspins ...
    3285             : !> \param xalmo_history ...
    3286             : !> \param assume_t0_q0x ...
    3287             : !> \param optimize_theta ...
    3288             : !> \param envelope_amplitude ...
    3289             : !> \param eps_filter ...
    3290             : !> \param order_lanczos ...
    3291             : !> \param eps_lanczos ...
    3292             : !> \param max_iter_lanczos ...
    3293             : !> \param nocc_of_domain ...
    3294             : !> \par History
    3295             : !>       2016.11 created [Rustam Z Khaliullin]
    3296             : !> \author Rustam Z Khaliullin
    3297             : ! **************************************************************************************************
    3298         104 :    SUBROUTINE xalmo_initial_guess(m_guess, m_t_in, m_t0, m_quench_t, &
    3299         104 :                                   m_overlap, m_sigma_tmpl, nspins, xalmo_history, assume_t0_q0x, &
    3300             :                                   optimize_theta, envelope_amplitude, eps_filter, order_lanczos, eps_lanczos, &
    3301         104 :                                   max_iter_lanczos, nocc_of_domain)
    3302             : 
    3303             :       TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT)      :: m_guess
    3304             :       TYPE(dbcsr_type), DIMENSION(:), INTENT(IN)         :: m_t_in, m_t0, m_quench_t
    3305             :       TYPE(dbcsr_type), INTENT(IN)                       :: m_overlap
    3306             :       TYPE(dbcsr_type), DIMENSION(:), INTENT(IN)         :: m_sigma_tmpl
    3307             :       INTEGER, INTENT(IN)                                :: nspins
    3308             :       TYPE(almo_scf_history_type), INTENT(IN)            :: xalmo_history
    3309             :       LOGICAL, INTENT(IN)                                :: assume_t0_q0x, optimize_theta
    3310             :       REAL(KIND=dp), INTENT(IN)                          :: envelope_amplitude, eps_filter
    3311             :       INTEGER, INTENT(IN)                                :: order_lanczos
    3312             :       REAL(KIND=dp), INTENT(IN)                          :: eps_lanczos
    3313             :       INTEGER, INTENT(IN)                                :: max_iter_lanczos
    3314             :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: nocc_of_domain
    3315             : 
    3316             :       CHARACTER(len=*), PARAMETER :: routineN = 'xalmo_initial_guess'
    3317             : 
    3318             :       INTEGER                                            :: handle, iaspc, ispin, istore, naspc, &
    3319             :                                                             unit_nr
    3320             :       LOGICAL                                            :: aspc_guess
    3321             :       REAL(KIND=dp)                                      :: alpha
    3322             :       TYPE(cp_logger_type), POINTER                      :: logger
    3323             :       TYPE(dbcsr_type)                                   :: m_extrapolated, m_sigma_tmp
    3324             : 
    3325         104 :       CALL timeset(routineN, handle)
    3326             : 
    3327             :       ! get a useful output_unit
    3328         104 :       logger => cp_get_default_logger()
    3329         104 :       IF (logger%para_env%is_source()) THEN
    3330          52 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    3331             :       ELSE
    3332             :          unit_nr = -1
    3333             :       END IF
    3334             : 
    3335         104 :       IF (optimize_theta) THEN
    3336           0 :          CPWARN("unused option")
    3337             :          ! just not to trigger unused variable
    3338           0 :          alpha = envelope_amplitude
    3339             :       END IF
    3340             : 
    3341             :       ! if extrapolation order is zero then the standard guess is used
    3342             :       ! ... the number of stored history points will remain zero if extrapolation order is zero
    3343         104 :       IF (xalmo_history%istore == 0) THEN
    3344             :          aspc_guess = .FALSE.
    3345             :       ELSE
    3346             :          aspc_guess = .TRUE.
    3347             :       END IF
    3348             : 
    3349             :       ! create initial guess
    3350             :       IF (.NOT. aspc_guess) THEN
    3351             : 
    3352         124 :          DO ispin = 1, nspins
    3353             : 
    3354             :             ! zero initial guess for the delocalization amplitudes
    3355             :             ! or the supplied guess for orbitals
    3356         124 :             IF (assume_t0_q0x) THEN
    3357          30 :                CALL dbcsr_set(m_guess(ispin), 0.0_dp)
    3358             :             ELSE
    3359             :                ! copy coefficients to m_guess
    3360          32 :                CALL dbcsr_copy(m_guess(ispin), m_t_in(ispin))
    3361             :             END IF
    3362             : 
    3363             :          END DO !ispins
    3364             : 
    3365             :       ELSE !aspc_guess
    3366             : 
    3367          42 :          CALL cite_reference(Kolafa2004)
    3368          42 :          CALL cite_reference(Kuhne2007)
    3369             : 
    3370          42 :          naspc = MIN(xalmo_history%istore, xalmo_history%nstore)
    3371          42 :          IF (unit_nr > 0) THEN
    3372             :             WRITE (unit_nr, FMT="(/,T2,A,/,/,T3,A,I0)") &
    3373          21 :                "Parameters for the always stable predictor-corrector (ASPC) method:", &
    3374          42 :                "ASPC order: ", naspc
    3375             :          END IF
    3376             : 
    3377          84 :          DO ispin = 1, nspins
    3378             : 
    3379             :             CALL dbcsr_create(m_extrapolated, &
    3380          42 :                               template=m_quench_t(ispin), matrix_type=dbcsr_type_no_symmetry)
    3381             :             CALL dbcsr_create(m_sigma_tmp, &
    3382          42 :                               template=m_sigma_tmpl(ispin), matrix_type=dbcsr_type_no_symmetry)
    3383             : 
    3384             :             ! set to zero before accumulation
    3385          42 :             CALL dbcsr_set(m_guess(ispin), 0.0_dp)
    3386             : 
    3387             :             ! extrapolation
    3388         124 :             DO iaspc = 1, naspc
    3389             : 
    3390          82 :                istore = MOD(xalmo_history%istore - iaspc, xalmo_history%nstore) + 1
    3391             :                alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
    3392          82 :                        binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
    3393          82 :                IF (unit_nr > 0) THEN
    3394             :                   WRITE (unit_nr, FMT="(T3,A2,I0,A4,F10.6)") &
    3395          41 :                      "B(", iaspc, ") = ", alpha
    3396             :                END IF
    3397             : 
    3398             :                ! m_extrapolated - initialize the correct sparsity pattern
    3399             :                ! it must be kept throughout extrapolation
    3400          82 :                CALL dbcsr_copy(m_extrapolated, m_quench_t(ispin))
    3401             : 
    3402             :                ! project t0 onto the previous DMs
    3403             :                ! note that t0 is projected instead of any other matrix (e.g.
    3404             :                ! t_SCF from the prev step or random t)
    3405             :                ! this is done to keep orbitals phase (i.e. sign) the same as in
    3406             :                ! t0. if this is not done then subtracting t0 on the next step
    3407             :                ! will produce a terrible guess and extrapolation will fail
    3408             :                CALL dbcsr_multiply("N", "N", 1.0_dp, &
    3409             :                                    xalmo_history%matrix_p_up_down(ispin, istore), &
    3410             :                                    m_t0(ispin), &
    3411             :                                    0.0_dp, m_extrapolated, &
    3412          82 :                                    retain_sparsity=.TRUE.)
    3413             :                ! normalize MOs
    3414             :                CALL orthogonalize_mos(ket=m_extrapolated, &
    3415             :                                       overlap=m_sigma_tmp, &
    3416             :                                       metric=m_overlap, &
    3417             :                                       retain_locality=.TRUE., &
    3418             :                                       only_normalize=.FALSE., &
    3419             :                                       nocc_of_domain=nocc_of_domain(:, ispin), &
    3420             :                                       eps_filter=eps_filter, &
    3421             :                                       order_lanczos=order_lanczos, &
    3422             :                                       eps_lanczos=eps_lanczos, &
    3423          82 :                                       max_iter_lanczos=max_iter_lanczos)
    3424             : 
    3425             :                ! now accumulate. correct sparsity is ensured
    3426             :                CALL dbcsr_add(m_guess(ispin), m_extrapolated, &
    3427         124 :                               1.0_dp, (1.0_dp*alpha)/naspc)
    3428             : 
    3429             :             END DO !iaspc
    3430             : 
    3431          42 :             CALL dbcsr_release(m_extrapolated)
    3432             : 
    3433             :             ! normalize MOs
    3434             :             CALL orthogonalize_mos(ket=m_guess(ispin), &
    3435             :                                    overlap=m_sigma_tmp, &
    3436             :                                    metric=m_overlap, &
    3437             :                                    retain_locality=.TRUE., &
    3438             :                                    only_normalize=.FALSE., &
    3439             :                                    nocc_of_domain=nocc_of_domain(:, ispin), &
    3440             :                                    eps_filter=eps_filter, &
    3441             :                                    order_lanczos=order_lanczos, &
    3442             :                                    eps_lanczos=eps_lanczos, &
    3443          42 :                                    max_iter_lanczos=max_iter_lanczos)
    3444             : 
    3445          42 :             CALL dbcsr_release(m_sigma_tmp)
    3446             : 
    3447             :             ! project the t0 space out from the extrapolated state
    3448             :             ! this can be done outside this subroutine
    3449          84 :             IF (assume_t0_q0x) THEN
    3450             :                CALL dbcsr_add(m_guess(ispin), m_t0(ispin), &
    3451          12 :                               1.0_dp, -1.0_dp)
    3452             :             END IF !assume_t0_q0x
    3453             : 
    3454             :          END DO !ispin
    3455             : 
    3456             :       END IF !aspc_guess?
    3457             : 
    3458         104 :       CALL timestop(handle)
    3459             : 
    3460         104 :    END SUBROUTINE xalmo_initial_guess
    3461             : 
    3462             : END MODULE almo_scf_methods
    3463             : 

Generated by: LCOV version 1.15