LCOV - code coverage report
Current view: top level - src - ct_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 0 563 0.0 %
Date: 2025-02-18 08:24:35 Functions: 0 9 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Cayley transformation methods
      10             : !> \par History
      11             : !>       2011.06 created [Rustam Z Khaliullin]
      12             : !> \author Rustam Z Khaliullin
      13             : ! **************************************************************************************************
      14             : MODULE ct_methods
      15             :    USE cp_dbcsr_api,                    ONLY: &
      16             :         dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, &
      17             :         dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, dbcsr_finalize, &
      18             :         dbcsr_get_diag, dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
      19             :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      20             :         dbcsr_multiply, dbcsr_nblkcols_total, dbcsr_nblkrows_total, dbcsr_release, &
      21             :         dbcsr_reserve_block2d, dbcsr_scale, dbcsr_set, dbcsr_set_diag, dbcsr_transposed, &
      22             :         dbcsr_type, dbcsr_type_no_symmetry, dbcsr_work_create
      23             :    USE cp_dbcsr_cholesky,               ONLY: cp_dbcsr_cholesky_decompose,&
      24             :                                               cp_dbcsr_cholesky_invert
      25             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_frobenius_norm,&
      26             :                                               dbcsr_hadamard_product,&
      27             :                                               dbcsr_maxabs
      28             :    USE cp_dbcsr_diag,                   ONLY: cp_dbcsr_syevd
      29             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      30             :                                               cp_logger_get_default_unit_nr,&
      31             :                                               cp_logger_type
      32             :    USE ct_types,                        ONLY: ct_step_env_type
      33             :    USE input_constants,                 ONLY: &
      34             :         cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, &
      35             :         cg_liu_storey, cg_polak_ribiere, cg_zero, tensor_orthogonal, tensor_up_down
      36             :    USE iterate_matrix,                  ONLY: matrix_sqrt_Newton_Schulz
      37             :    USE kinds,                           ONLY: dp
      38             :    USE machine,                         ONLY: m_walltime
      39             : #include "./base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ct_methods'
      46             : 
      47             :    ! Public subroutines
      48             :    PUBLIC :: ct_step_execute, analytic_line_search, diagonalize_diagonal_blocks
      49             : 
      50             : CONTAINS
      51             : 
      52             : ! **************************************************************************************************
      53             : !> \brief Performs Cayley transformation
      54             : !> \param cts_env ...
      55             : !> \par History
      56             : !>       2011.06 created [Rustam Z Khaliullin]
      57             : !> \author Rustam Z Khaliullin
      58             : ! **************************************************************************************************
      59           0 :    SUBROUTINE ct_step_execute(cts_env)
      60             : 
      61             :       TYPE(ct_step_env_type)                             :: cts_env
      62             : 
      63             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ct_step_execute'
      64             : 
      65             :       INTEGER                                            :: handle, n, preconditioner_type, unit_nr
      66             :       REAL(KIND=dp)                                      :: gap_estimate, safety_margin
      67           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
      68             :       TYPE(cp_logger_type), POINTER                      :: logger
      69             :       TYPE(dbcsr_type)                                   :: matrix_pp, matrix_pq, matrix_qp, &
      70             :                                                             matrix_qp_save, matrix_qq, oo1, &
      71             :                                                             oo1_sqrt, oo1_sqrt_inv, t_corr, tmp1, &
      72             :                                                             u_pp, u_qq
      73             : 
      74             : !TYPE(dbcsr_type)                :: rst_x1, rst_x2
      75             : !REAL(KIND=dp)                      :: ener_tmp
      76             : !TYPE(dbcsr_iterator_type)            :: iter
      77             : !INTEGER                            :: iblock_row,iblock_col,&
      78             : !                                      iblock_row_size,iblock_col_size
      79             : !REAL(KIND=dp), DIMENSION(:,:), POINTER :: data_p
      80             : 
      81           0 :       CALL timeset(routineN, handle)
      82             : 
      83           0 :       logger => cp_get_default_logger()
      84           0 :       IF (logger%para_env%is_source()) THEN
      85           0 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
      86             :       ELSE
      87             :          unit_nr = -1
      88             :       END IF
      89             : 
      90             :       ! check if all input is in place and flags are consistent
      91           0 :       IF (cts_env%update_q .AND. (.NOT. cts_env%update_p)) THEN
      92           0 :          CPABORT("q-update is possible only with p-update")
      93             :       END IF
      94             : 
      95           0 :       IF (cts_env%tensor_type .EQ. tensor_up_down) THEN
      96           0 :          CPABORT("riccati is not implemented for biorthogonal basis")
      97             :       END IF
      98             : 
      99           0 :       IF (.NOT. ASSOCIATED(cts_env%matrix_ks)) THEN
     100           0 :          CPABORT("KS matrix is not associated")
     101             :       END IF
     102             : 
     103           0 :       IF (cts_env%use_virt_orbs .AND. (.NOT. cts_env%use_occ_orbs)) THEN
     104           0 :          CPABORT("virtual orbs can be used only with occupied orbs")
     105             :       END IF
     106             : 
     107           0 :       IF (cts_env%use_occ_orbs) THEN
     108           0 :          IF (.NOT. ASSOCIATED(cts_env%matrix_t)) THEN
     109           0 :             CPABORT("T matrix is not associated")
     110             :          END IF
     111           0 :          IF (.NOT. ASSOCIATED(cts_env%matrix_qp_template)) THEN
     112           0 :             CPABORT("QP template is not associated")
     113             :          END IF
     114           0 :          IF (.NOT. ASSOCIATED(cts_env%matrix_pq_template)) THEN
     115           0 :             CPABORT("PQ template is not associated")
     116             :          END IF
     117             :       END IF
     118             : 
     119           0 :       IF (cts_env%use_virt_orbs) THEN
     120           0 :          IF (.NOT. ASSOCIATED(cts_env%matrix_v)) THEN
     121           0 :             CPABORT("V matrix is not associated")
     122             :          END IF
     123             :       ELSE
     124           0 :          IF (.NOT. ASSOCIATED(cts_env%matrix_p)) THEN
     125           0 :             CPABORT("P matrix is not associated")
     126             :          END IF
     127             :       END IF
     128             : 
     129           0 :       IF (cts_env%tensor_type .NE. tensor_up_down .AND. &
     130             :           cts_env%tensor_type .NE. tensor_orthogonal) THEN
     131           0 :          CPABORT("illegal tensor flag")
     132             :       END IF
     133             : 
     134             :       ! start real calculations
     135           0 :       IF (cts_env%use_occ_orbs) THEN
     136             : 
     137             :          ! create matrices for various ks blocks
     138             :          CALL dbcsr_create(matrix_pp, &
     139             :                            template=cts_env%p_index_up, &
     140           0 :                            matrix_type=dbcsr_type_no_symmetry)
     141             :          CALL dbcsr_create(matrix_qp, &
     142             :                            template=cts_env%matrix_qp_template, &
     143           0 :                            matrix_type=dbcsr_type_no_symmetry)
     144             :          CALL dbcsr_create(matrix_qq, &
     145             :                            template=cts_env%q_index_up, &
     146           0 :                            matrix_type=dbcsr_type_no_symmetry)
     147             :          CALL dbcsr_create(matrix_pq, &
     148             :                            template=cts_env%matrix_pq_template, &
     149           0 :                            matrix_type=dbcsr_type_no_symmetry)
     150             : 
     151             :          ! create the residue matrix
     152             :          CALL dbcsr_create(cts_env%matrix_res, &
     153           0 :                            template=cts_env%matrix_qp_template)
     154             : 
     155             :          CALL assemble_ks_qp_blocks(cts_env%matrix_ks, &
     156             :                                     cts_env%matrix_p, &
     157             :                                     cts_env%matrix_t, &
     158             :                                     cts_env%matrix_v, &
     159             :                                     cts_env%q_index_down, &
     160             :                                     cts_env%p_index_up, &
     161             :                                     cts_env%q_index_up, &
     162             :                                     matrix_pp, &
     163             :                                     matrix_qq, &
     164             :                                     matrix_qp, &
     165             :                                     matrix_pq, &
     166             :                                     cts_env%tensor_type, &
     167             :                                     cts_env%use_virt_orbs, &
     168           0 :                                     cts_env%eps_filter)
     169             : 
     170             :          ! create a matrix of single-excitation amplitudes
     171             :          CALL dbcsr_create(cts_env%matrix_x, &
     172           0 :                            template=cts_env%matrix_qp_template)
     173           0 :          IF (ASSOCIATED(cts_env%matrix_x_guess)) THEN
     174             :             CALL dbcsr_copy(cts_env%matrix_x, &
     175           0 :                             cts_env%matrix_x_guess)
     176           0 :             IF (cts_env%tensor_type .EQ. tensor_orthogonal) THEN
     177             :                ! bring x from contravariant-covariant representation
     178             :                ! to the orthogonal/cholesky representation
     179             :                ! use res as temporary storage
     180             :                CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_down, &
     181             :                                    cts_env%matrix_x, 0.0_dp, cts_env%matrix_res, &
     182           0 :                                    filter_eps=cts_env%eps_filter)
     183             :                CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_res, &
     184             :                                    cts_env%p_index_up, 0.0_dp, &
     185             :                                    cts_env%matrix_x, &
     186           0 :                                    filter_eps=cts_env%eps_filter)
     187             :             END IF
     188             :          ELSE
     189             :             ! set amplitudes to zero
     190           0 :             CALL dbcsr_set(cts_env%matrix_x, 0.0_dp)
     191             :          END IF
     192             : 
     193             :          !SELECT CASE (cts_env%preconditioner_type)
     194             :          !CASE (prec_eigenvector_blocks,prec_eigenvector_full)
     195           0 :          preconditioner_type = 1
     196           0 :          safety_margin = 2.0_dp
     197           0 :          gap_estimate = 0.0001_dp
     198             :          SELECT CASE (preconditioner_type)
     199             :          CASE (1, 2)
     200             : !RZK-warning diagonalization works only with orthogonal tensor!!!
     201             :             ! find a better basis by diagonalizing diagonal blocks
     202             :             ! first pp
     203             :             CALL dbcsr_create(u_pp, template=matrix_pp, &
     204           0 :                               matrix_type=dbcsr_type_no_symmetry)
     205             :             !IF (cts_env%preconditioner_type.eq.prec_eigenvector_full) THEN
     206             :             IF (.TRUE.) THEN
     207           0 :                CALL dbcsr_get_info(matrix_pp, nfullrows_total=n)
     208           0 :                ALLOCATE (evals(n))
     209             :                CALL cp_dbcsr_syevd(matrix_pp, u_pp, evals, &
     210           0 :                                    cts_env%para_env, cts_env%blacs_env)
     211           0 :                DEALLOCATE (evals)
     212             :             ELSE
     213             :                CALL diagonalize_diagonal_blocks(matrix_pp, u_pp)
     214             :             END IF
     215             :             ! and now qq
     216             :             CALL dbcsr_create(u_qq, template=matrix_qq, &
     217           0 :                               matrix_type=dbcsr_type_no_symmetry)
     218             :             !IF (cts_env%preconditioner_type.eq.prec_eigenvector_full) THEN
     219             :             IF (.TRUE.) THEN
     220           0 :                CALL dbcsr_get_info(matrix_qq, nfullrows_total=n)
     221           0 :                ALLOCATE (evals(n))
     222             :                CALL cp_dbcsr_syevd(matrix_qq, u_qq, evals, &
     223           0 :                                    cts_env%para_env, cts_env%blacs_env)
     224           0 :                DEALLOCATE (evals)
     225             :             ELSE
     226             :                CALL diagonalize_diagonal_blocks(matrix_qq, u_qq)
     227             :             END IF
     228             : 
     229             :             ! apply the transformation to all matrices
     230             :             CALL matrix_forward_transform(matrix_pp, u_pp, u_pp, &
     231           0 :                                           cts_env%eps_filter)
     232             :             CALL matrix_forward_transform(matrix_qq, u_qq, u_qq, &
     233           0 :                                           cts_env%eps_filter)
     234             :             CALL matrix_forward_transform(matrix_qp, u_qq, u_pp, &
     235           0 :                                           cts_env%eps_filter)
     236             :             CALL matrix_forward_transform(matrix_pq, u_pp, u_qq, &
     237           0 :                                           cts_env%eps_filter)
     238             :             CALL matrix_forward_transform(cts_env%matrix_x, u_qq, u_pp, &
     239           0 :                                           cts_env%eps_filter)
     240             : 
     241           0 :             IF (cts_env%max_iter .GE. 0) THEN
     242             : 
     243             :                CALL solve_riccati_equation( &
     244             :                   pp=matrix_pp, &
     245             :                   qq=matrix_qq, &
     246             :                   qp=matrix_qp, &
     247             :                   pq=matrix_pq, &
     248             :                   x=cts_env%matrix_x, &
     249             :                   res=cts_env%matrix_res, &
     250             :                   neglect_quadratic_term=cts_env%neglect_quadratic_term, &
     251             :                   conjugator=cts_env%conjugator, &
     252             :                   max_iter=cts_env%max_iter, &
     253             :                   eps_convergence=cts_env%eps_convergence, &
     254             :                   eps_filter=cts_env%eps_filter, &
     255           0 :                   converged=cts_env%converged)
     256             : 
     257           0 :                IF (cts_env%converged) THEN
     258             :                   !IF (unit_nr>0) THEN
     259             :                   !   WRITE(unit_nr,*)
     260             :                   !   WRITE(unit_nr,'(T6,A)') &
     261             :                   !         "RICCATI equations solved"
     262             :                   !   CALL m_flush(unit_nr)
     263             :                   !ENDIF
     264             :                ELSE
     265           0 :                   CPABORT("RICCATI: CG algorithm has NOT converged")
     266             :                END IF
     267             : 
     268             :             END IF
     269             : 
     270           0 :             IF (cts_env%calculate_energy_corr) THEN
     271             : 
     272           0 :                CALL dbcsr_dot(matrix_qp, cts_env%matrix_x, cts_env%energy_correction)
     273             : 
     274             :             END IF
     275             : 
     276           0 :             CALL dbcsr_release(matrix_pp)
     277           0 :             CALL dbcsr_release(matrix_qp)
     278           0 :             CALL dbcsr_release(matrix_qq)
     279           0 :             CALL dbcsr_release(matrix_pq)
     280             : 
     281             :             ! back-transform to the original basis
     282             :             CALL matrix_backward_transform(cts_env%matrix_x, u_qq, &
     283           0 :                                            u_pp, cts_env%eps_filter)
     284             : 
     285           0 :             CALL dbcsr_release(u_qq)
     286           0 :             CALL dbcsr_release(u_pp)
     287             : 
     288             :             !CASE (prec_cholesky_inverse)
     289             :          CASE (3)
     290             : 
     291             : ! RZK-warning implemented only for orthogonal tensors!!!
     292             : ! generalization to up_down should be easy
     293             :             CALL dbcsr_create(u_pp, template=matrix_pp, &
     294             :                               matrix_type=dbcsr_type_no_symmetry)
     295             :             CALL dbcsr_copy(u_pp, matrix_pp)
     296             :             CALL dbcsr_scale(u_pp, -1.0_dp)
     297             :             CALL dbcsr_add_on_diag(u_pp, &
     298             :                                    ABS(safety_margin*gap_estimate))
     299             :             CALL cp_dbcsr_cholesky_decompose(u_pp, &
     300             :                                              para_env=cts_env%para_env, &
     301             :                                              blacs_env=cts_env%blacs_env)
     302             :             CALL cp_dbcsr_cholesky_invert(u_pp, &
     303             :                                           para_env=cts_env%para_env, &
     304             :                                           blacs_env=cts_env%blacs_env, &
     305             :                                           uplo_to_full=.TRUE.)
     306             :             !CALL dbcsr_scale(u_pp,-1.0_dp)
     307             : 
     308             :             CALL dbcsr_create(u_qq, template=matrix_qq, &
     309             :                               matrix_type=dbcsr_type_no_symmetry)
     310             :             CALL dbcsr_copy(u_qq, matrix_qq)
     311             :             CALL dbcsr_add_on_diag(u_qq, &
     312             :                                    ABS(safety_margin*gap_estimate))
     313             :             CALL cp_dbcsr_cholesky_decompose(u_qq, &
     314             :                                              para_env=cts_env%para_env, &
     315             :                                              blacs_env=cts_env%blacs_env)
     316             :             CALL cp_dbcsr_cholesky_invert(u_qq, &
     317             :                                           para_env=cts_env%para_env, &
     318             :                                           blacs_env=cts_env%blacs_env, &
     319             :                                           uplo_to_full=.TRUE.)
     320             : 
     321             :             ! transform all riccati matrices (left-right preconditioner)
     322             :             CALL dbcsr_create(tmp1, template=matrix_qq, &
     323             :                               matrix_type=dbcsr_type_no_symmetry)
     324             :             CALL dbcsr_multiply("N", "N", 1.0_dp, u_qq, &
     325             :                                 matrix_qq, 0.0_dp, tmp1, &
     326             :                                 filter_eps=cts_env%eps_filter)
     327             :             CALL dbcsr_copy(matrix_qq, tmp1)
     328             :             CALL dbcsr_release(tmp1)
     329             : 
     330             :             CALL dbcsr_create(tmp1, template=matrix_pp, &
     331             :                               matrix_type=dbcsr_type_no_symmetry)
     332             :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_pp, &
     333             :                                 u_pp, 0.0_dp, tmp1, &
     334             :                                 filter_eps=cts_env%eps_filter)
     335             :             CALL dbcsr_copy(matrix_pp, tmp1)
     336             :             CALL dbcsr_release(tmp1)
     337             : 
     338             :             CALL dbcsr_create(matrix_qp_save, template=matrix_qp, &
     339             :                               matrix_type=dbcsr_type_no_symmetry)
     340             :             CALL dbcsr_copy(matrix_qp_save, matrix_qp)
     341             : 
     342             :             CALL dbcsr_create(tmp1, template=matrix_qp, &
     343             :                               matrix_type=dbcsr_type_no_symmetry)
     344             :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
     345             :                                 u_pp, 0.0_dp, tmp1, &
     346             :                                 filter_eps=cts_env%eps_filter)
     347             :             CALL dbcsr_multiply("N", "N", 1.0_dp, u_qq, tmp1, &
     348             :                                 0.0_dp, matrix_qp, &
     349             :                                 filter_eps=cts_env%eps_filter)
     350             :             CALL dbcsr_release(tmp1)
     351             : !CALL dbcsr_print(matrix_qq)
     352             : !CALL dbcsr_print(matrix_qp)
     353             : !CALL dbcsr_print(matrix_pp)
     354             : 
     355             :             IF (cts_env%max_iter .GE. 0) THEN
     356             : 
     357             :                CALL solve_riccati_equation( &
     358             :                   pp=matrix_pp, &
     359             :                   qq=matrix_qq, &
     360             :                   qp=matrix_qp, &
     361             :                   pq=matrix_pq, &
     362             :                   oo=u_pp, &
     363             :                   vv=u_qq, &
     364             :                   x=cts_env%matrix_x, &
     365             :                   res=cts_env%matrix_res, &
     366             :                   neglect_quadratic_term=cts_env%neglect_quadratic_term, &
     367             :                   conjugator=cts_env%conjugator, &
     368             :                   max_iter=cts_env%max_iter, &
     369             :                   eps_convergence=cts_env%eps_convergence, &
     370             :                   eps_filter=cts_env%eps_filter, &
     371             :                   converged=cts_env%converged)
     372             : 
     373             :                IF (cts_env%converged) THEN
     374             :                   !IF (unit_nr>0) THEN
     375             :                   !   WRITE(unit_nr,*)
     376             :                   !   WRITE(unit_nr,'(T6,A)') &
     377             :                   !         "RICCATI equations solved"
     378             :                   !   CALL m_flush(unit_nr)
     379             :                   !ENDIF
     380             :                ELSE
     381             :                   CPABORT("RICCATI: CG algorithm has NOT converged")
     382             :                END IF
     383             : 
     384             :             END IF
     385             : 
     386             :             IF (cts_env%calculate_energy_corr) THEN
     387             : 
     388             :                CALL dbcsr_dot(matrix_qp_save, cts_env%matrix_x, cts_env%energy_correction)
     389             : 
     390             :             END IF
     391             :             CALL dbcsr_release(matrix_qp_save)
     392             : 
     393             :             CALL dbcsr_release(matrix_pp)
     394             :             CALL dbcsr_release(matrix_qp)
     395             :             CALL dbcsr_release(matrix_qq)
     396             :             CALL dbcsr_release(matrix_pq)
     397             : 
     398             :             CALL dbcsr_release(u_qq)
     399             :             CALL dbcsr_release(u_pp)
     400             : 
     401             :          CASE DEFAULT
     402             :             CPABORT("illegal preconditioner type")
     403             :          END SELECT ! preconditioner type
     404             : 
     405           0 :          IF (cts_env%update_p) THEN
     406             : 
     407           0 :             IF (cts_env%tensor_type .EQ. tensor_up_down) THEN
     408           0 :                CPABORT("orbital update is NYI for this tensor type")
     409             :             END IF
     410             : 
     411             :             ! transform occupied orbitals
     412             :             ! in a way that preserves the overlap metric
     413             :             CALL dbcsr_create(oo1, &
     414             :                               template=cts_env%p_index_up, &
     415           0 :                               matrix_type=dbcsr_type_no_symmetry)
     416             :             CALL dbcsr_create(oo1_sqrt_inv, &
     417           0 :                               template=oo1)
     418             :             CALL dbcsr_create(oo1_sqrt, &
     419           0 :                               template=oo1)
     420             : 
     421             :             ! Compute (1+tr(X).X)^(-1/2)_up_down
     422             :             CALL dbcsr_multiply("T", "N", 1.0_dp, cts_env%matrix_x, &
     423             :                                 cts_env%matrix_x, 0.0_dp, oo1, &
     424           0 :                                 filter_eps=cts_env%eps_filter)
     425           0 :             CALL dbcsr_add_on_diag(oo1, 1.0_dp)
     426             :             CALL matrix_sqrt_Newton_Schulz(oo1_sqrt, &
     427             :                                            oo1_sqrt_inv, &
     428             :                                            oo1, &
     429             :                                            !if cholesky is used then sqrt
     430             :                                            !guess cannot be provided
     431             :                                            !matrix_sqrt_inv_guess=cts_env%p_index_up,&
     432             :                                            !matrix_sqrt_guess=cts_env%p_index_down,&
     433             :                                            threshold=cts_env%eps_filter, &
     434             :                                            order=cts_env%order_lanczos, &
     435             :                                            eps_lanczos=cts_env%eps_lancsoz, &
     436           0 :                                            max_iter_lanczos=cts_env%max_iter_lanczos)
     437             :             CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%p_index_up, &
     438             :                                 oo1_sqrt_inv, 0.0_dp, oo1, &
     439           0 :                                 filter_eps=cts_env%eps_filter)
     440             :             CALL dbcsr_multiply("N", "N", 1.0_dp, oo1, &
     441             :                                 cts_env%p_index_down, 0.0_dp, oo1_sqrt, &
     442           0 :                                 filter_eps=cts_env%eps_filter)
     443           0 :             CALL dbcsr_release(oo1)
     444           0 :             CALL dbcsr_release(oo1_sqrt_inv)
     445             : 
     446             :             ! bring x to contravariant-covariant representation now
     447             :             CALL dbcsr_create(matrix_qp, &
     448             :                               template=cts_env%matrix_qp_template, &
     449           0 :                               matrix_type=dbcsr_type_no_symmetry)
     450             :             CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_up, &
     451             :                                 cts_env%matrix_x, 0.0_dp, matrix_qp, &
     452           0 :                                 filter_eps=cts_env%eps_filter)
     453             :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
     454             :                                 cts_env%p_index_down, 0.0_dp, &
     455             :                                 cts_env%matrix_x, &
     456           0 :                                 filter_eps=cts_env%eps_filter)
     457           0 :             CALL dbcsr_release(matrix_qp)
     458             : 
     459             :             ! update T=T+X or T=T+V.X (whichever is appropriate)
     460           0 :             CALL dbcsr_create(t_corr, template=cts_env%matrix_t)
     461           0 :             IF (cts_env%use_virt_orbs) THEN
     462             :                CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_v, &
     463             :                                    cts_env%matrix_x, 0.0_dp, t_corr, &
     464           0 :                                    filter_eps=cts_env%eps_filter)
     465             :                CALL dbcsr_add(cts_env%matrix_t, t_corr, &
     466           0 :                               1.0_dp, 1.0_dp)
     467             :             ELSE
     468             :                CALL dbcsr_add(cts_env%matrix_t, cts_env%matrix_x, &
     469           0 :                               1.0_dp, 1.0_dp)
     470             :             END IF
     471             :             ! adjust T so the metric is preserved: T=(T+X).(1+tr(X).X)^(-1/2)
     472             :             CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_t, oo1_sqrt, &
     473           0 :                                 0.0_dp, t_corr, filter_eps=cts_env%eps_filter)
     474           0 :             CALL dbcsr_copy(cts_env%matrix_t, t_corr)
     475             : 
     476           0 :             CALL dbcsr_release(t_corr)
     477           0 :             CALL dbcsr_release(oo1_sqrt)
     478             : 
     479             :          ELSE ! do not update p
     480             : 
     481           0 :             IF (cts_env%tensor_type .EQ. tensor_orthogonal) THEN
     482             :                ! bring x to contravariant-covariant representation
     483             :                CALL dbcsr_create(matrix_qp, &
     484             :                                  template=cts_env%matrix_qp_template, &
     485           0 :                                  matrix_type=dbcsr_type_no_symmetry)
     486             :                CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_up, &
     487             :                                    cts_env%matrix_x, 0.0_dp, matrix_qp, &
     488           0 :                                    filter_eps=cts_env%eps_filter)
     489             :                CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
     490             :                                    cts_env%p_index_down, 0.0_dp, &
     491             :                                    cts_env%matrix_x, &
     492           0 :                                    filter_eps=cts_env%eps_filter)
     493           0 :                CALL dbcsr_release(matrix_qp)
     494             :             END IF
     495             : 
     496             :          END IF
     497             : 
     498             :       ELSE
     499           0 :          CPABORT("illegal occ option")
     500             :       END IF
     501             : 
     502           0 :       CALL timestop(handle)
     503             : 
     504           0 :    END SUBROUTINE ct_step_execute
     505             : 
     506             : ! **************************************************************************************************
     507             : !> \brief computes oo, ov, vo, and vv blocks of the ks matrix
     508             : !> \param ks ...
     509             : !> \param p ...
     510             : !> \param t ...
     511             : !> \param v ...
     512             : !> \param q_index_down ...
     513             : !> \param p_index_up ...
     514             : !> \param q_index_up ...
     515             : !> \param pp ...
     516             : !> \param qq ...
     517             : !> \param qp ...
     518             : !> \param pq ...
     519             : !> \param tensor_type ...
     520             : !> \param use_virt_orbs ...
     521             : !> \param eps_filter ...
     522             : !> \par History
     523             : !>       2011.06 created [Rustam Z Khaliullin]
     524             : !> \author Rustam Z Khaliullin
     525             : ! **************************************************************************************************
     526           0 :    SUBROUTINE assemble_ks_qp_blocks(ks, p, t, v, q_index_down, &
     527             :                                     p_index_up, q_index_up, pp, qq, qp, pq, tensor_type, use_virt_orbs, eps_filter)
     528             : 
     529             :       TYPE(dbcsr_type), INTENT(IN)                       :: ks, p, t, v, q_index_down, p_index_up, &
     530             :                                                             q_index_up
     531             :       TYPE(dbcsr_type), INTENT(OUT)                      :: pp, qq, qp, pq
     532             :       INTEGER, INTENT(IN)                                :: tensor_type
     533             :       LOGICAL, INTENT(IN)                                :: use_virt_orbs
     534             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
     535             : 
     536             :       CHARACTER(len=*), PARAMETER :: routineN = 'assemble_ks_qp_blocks'
     537             : 
     538             :       INTEGER                                            :: handle
     539             :       LOGICAL                                            :: library_fixed
     540             :       TYPE(dbcsr_type)                                   :: kst, ksv, no, on, oo, q_index_up_nosym, &
     541             :                                                             sp, spf, t_or, v_or
     542             : 
     543           0 :       CALL timeset(routineN, handle)
     544             : 
     545           0 :       IF (use_virt_orbs) THEN
     546             : 
     547             :          ! orthogonalize the orbitals
     548           0 :          CALL dbcsr_create(t_or, template=t)
     549           0 :          CALL dbcsr_create(v_or, template=v)
     550             :          CALL dbcsr_multiply("N", "N", 1.0_dp, t, p_index_up, &
     551           0 :                              0.0_dp, t_or, filter_eps=eps_filter)
     552             :          CALL dbcsr_multiply("N", "N", 1.0_dp, v, q_index_up, &
     553           0 :                              0.0_dp, v_or, filter_eps=eps_filter)
     554             : 
     555             :          ! KS.T
     556           0 :          CALL dbcsr_create(kst, template=t)
     557             :          CALL dbcsr_multiply("N", "N", 1.0_dp, ks, t_or, &
     558           0 :                              0.0_dp, kst, filter_eps=eps_filter)
     559             :          ! pp=tr(T)*KS.T
     560             :          CALL dbcsr_multiply("T", "N", 1.0_dp, t_or, kst, &
     561           0 :                              0.0_dp, pp, filter_eps=eps_filter)
     562             :          ! qp=tr(V)*KS.T
     563             :          CALL dbcsr_multiply("T", "N", 1.0_dp, v_or, kst, &
     564           0 :                              0.0_dp, qp, filter_eps=eps_filter)
     565           0 :          CALL dbcsr_release(kst)
     566             : 
     567             :          ! KS.V
     568           0 :          CALL dbcsr_create(ksv, template=v)
     569             :          CALL dbcsr_multiply("N", "N", 1.0_dp, ks, v_or, &
     570           0 :                              0.0_dp, ksv, filter_eps=eps_filter)
     571             :          ! tr(T)*KS.V
     572             :          CALL dbcsr_multiply("T", "N", 1.0_dp, t_or, ksv, &
     573           0 :                              0.0_dp, pq, filter_eps=eps_filter)
     574             :          ! tr(V)*KS.V
     575             :          CALL dbcsr_multiply("T", "N", 1.0_dp, v_or, ksv, &
     576           0 :                              0.0_dp, qq, filter_eps=eps_filter)
     577           0 :          CALL dbcsr_release(ksv)
     578             : 
     579           0 :          CALL dbcsr_release(t_or)
     580           0 :          CALL dbcsr_release(v_or)
     581             : 
     582             :       ELSE ! no virtuals, use projected AOs
     583             : 
     584             : ! THIS PROCEDURE HAS NOT BEEN UPDATED FOR CHOLESKY p/q_index_up/down
     585             :          CALL dbcsr_create(sp, template=q_index_down, &
     586           0 :                            matrix_type=dbcsr_type_no_symmetry)
     587             :          CALL dbcsr_create(spf, template=q_index_down, &
     588           0 :                            matrix_type=dbcsr_type_no_symmetry)
     589             : 
     590             :          ! qp=KS*T
     591             :          CALL dbcsr_multiply("N", "N", 1.0_dp, ks, t, 0.0_dp, qp, &
     592           0 :                              filter_eps=eps_filter)
     593             :          ! pp=tr(T)*KS.T
     594             :          CALL dbcsr_multiply("T", "N", 1.0_dp, t, qp, 0.0_dp, pp, &
     595           0 :                              filter_eps=eps_filter)
     596             :          ! sp=-S_*P
     597             :          CALL dbcsr_multiply("N", "N", -1.0_dp, q_index_down, p, 0.0_dp, sp, &
     598           0 :                              filter_eps=eps_filter)
     599             : 
     600             :          ! sp=1/S^-S_.P
     601           0 :          SELECT CASE (tensor_type)
     602             :          CASE (tensor_up_down)
     603           0 :             CALL dbcsr_add_on_diag(sp, 1.0_dp)
     604             :          CASE (tensor_orthogonal)
     605             :             CALL dbcsr_create(q_index_up_nosym, template=q_index_up, &
     606           0 :                               matrix_type=dbcsr_type_no_symmetry)
     607           0 :             CALL dbcsr_desymmetrize(q_index_up, q_index_up_nosym)
     608           0 :             CALL dbcsr_add(sp, q_index_up_nosym, 1.0_dp, 1.0_dp)
     609           0 :             CALL dbcsr_release(q_index_up_nosym)
     610             :          END SELECT
     611             : 
     612             :          ! spf=(1/S^-S_.P)*KS
     613             :          CALL dbcsr_multiply("N", "N", 1.0_dp, sp, ks, 0.0_dp, spf, &
     614           0 :                              filter_eps=eps_filter)
     615             : 
     616             :          ! qp=spf*T
     617             :          CALL dbcsr_multiply("N", "N", 1.0_dp, spf, t, 0.0_dp, qp, &
     618           0 :                              filter_eps=eps_filter)
     619             : 
     620           0 :          SELECT CASE (tensor_type)
     621             :          CASE (tensor_up_down)
     622             :             ! pq=tr(qp)
     623           0 :             CALL dbcsr_transposed(pq, qp, transpose_distribution=.FALSE.)
     624             :          CASE (tensor_orthogonal)
     625             :             ! pq=sig^.tr(qp)
     626             :             CALL dbcsr_multiply("N", "T", 1.0_dp, p_index_up, qp, 0.0_dp, pq, &
     627           0 :                                 filter_eps=eps_filter)
     628           0 :             library_fixed = .FALSE.
     629           0 :             IF (library_fixed) THEN
     630             :                CALL dbcsr_transposed(qp, pq, transpose_distribution=.FALSE.)
     631             :             ELSE
     632             :                CALL dbcsr_create(no, template=qp, &
     633           0 :                                  matrix_type=dbcsr_type_no_symmetry)
     634             :                CALL dbcsr_multiply("N", "N", 1.0_dp, qp, p_index_up, 0.0_dp, no, &
     635           0 :                                    filter_eps=eps_filter)
     636           0 :                CALL dbcsr_copy(qp, no)
     637           0 :                CALL dbcsr_release(no)
     638             :             END IF
     639             :          END SELECT
     640             : 
     641             :          ! qq=spf*tr(sp)
     642             :          CALL dbcsr_multiply("N", "T", 1.0_dp, spf, sp, 0.0_dp, qq, &
     643           0 :                              filter_eps=eps_filter)
     644             : 
     645           0 :          SELECT CASE (tensor_type)
     646             :          CASE (tensor_up_down)
     647             : 
     648             :             CALL dbcsr_create(oo, template=pp, &
     649           0 :                               matrix_type=dbcsr_type_no_symmetry)
     650             :             CALL dbcsr_create(no, template=qp, &
     651           0 :                               matrix_type=dbcsr_type_no_symmetry)
     652             : 
     653             :             ! first index up
     654             :             CALL dbcsr_multiply("N", "N", 1.0_dp, q_index_up, qq, 0.0_dp, spf, &
     655           0 :                                 filter_eps=eps_filter)
     656           0 :             CALL dbcsr_copy(qq, spf)
     657             :             CALL dbcsr_multiply("N", "N", 1.0_dp, q_index_up, qp, 0.0_dp, no, &
     658           0 :                                 filter_eps=eps_filter)
     659           0 :             CALL dbcsr_copy(qp, no)
     660             :             CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pp, 0.0_dp, oo, &
     661           0 :                                 filter_eps=eps_filter)
     662           0 :             CALL dbcsr_copy(pp, oo)
     663             :             CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pq, 0.0_dp, on, &
     664           0 :                                 filter_eps=eps_filter)
     665           0 :             CALL dbcsr_copy(pq, on)
     666             : 
     667           0 :             CALL dbcsr_release(no)
     668           0 :             CALL dbcsr_release(oo)
     669             : 
     670             :          CASE (tensor_orthogonal)
     671             : 
     672             :             CALL dbcsr_create(oo, template=pp, &
     673           0 :                               matrix_type=dbcsr_type_no_symmetry)
     674             : 
     675             :             ! both indeces up in the pp block
     676             :             CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pp, 0.0_dp, oo, &
     677           0 :                                 filter_eps=eps_filter)
     678             :             CALL dbcsr_multiply("N", "N", 1.0_dp, oo, p_index_up, 0.0_dp, pp, &
     679           0 :                                 filter_eps=eps_filter)
     680             : 
     681           0 :             CALL dbcsr_release(oo)
     682             : 
     683             :          END SELECT
     684             : 
     685           0 :          CALL dbcsr_release(sp)
     686           0 :          CALL dbcsr_release(spf)
     687             : 
     688             :       END IF
     689             : 
     690           0 :       CALL timestop(handle)
     691             : 
     692           0 :    END SUBROUTINE assemble_ks_qp_blocks
     693             : 
     694             : ! **************************************************************************************************
     695             : !> \brief Solves the generalized Riccati or Sylvester eqation
     696             : !>        using the preconditioned conjugate gradient algorithm
     697             : !>          qp + qq.x.oo - vv.x.pp - vv.x.pq.x.oo = 0 [oo and vv are optional]
     698             : !>          qp + qq.x - x.pp - x.pq.x = 0
     699             : !> \param pp ...
     700             : !> \param qq ...
     701             : !> \param qp ...
     702             : !> \param pq ...
     703             : !> \param oo ...
     704             : !> \param vv ...
     705             : !> \param x ...
     706             : !> \param res ...
     707             : !> \param neglect_quadratic_term ...
     708             : !> \param conjugator ...
     709             : !> \param max_iter ...
     710             : !> \param eps_convergence ...
     711             : !> \param eps_filter ...
     712             : !> \param converged ...
     713             : !> \par History
     714             : !>       2011.06 created [Rustam Z Khaliullin]
     715             : !>       2011.11 generalized [Rustam Z Khaliullin]
     716             : !> \author Rustam Z Khaliullin
     717             : ! **************************************************************************************************
     718           0 :    RECURSIVE SUBROUTINE solve_riccati_equation(pp, qq, qp, pq, oo, vv, x, res, &
     719             :                                                neglect_quadratic_term, &
     720             :                                                conjugator, max_iter, eps_convergence, eps_filter, &
     721             :                                                converged)
     722             : 
     723             :       TYPE(dbcsr_type), INTENT(IN)                       :: pp, qq
     724             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: qp
     725             :       TYPE(dbcsr_type), INTENT(IN)                       :: pq
     726             :       TYPE(dbcsr_type), INTENT(IN), OPTIONAL             :: oo, vv
     727             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: x
     728             :       TYPE(dbcsr_type), INTENT(OUT)                      :: res
     729             :       LOGICAL, INTENT(IN)                                :: neglect_quadratic_term
     730             :       INTEGER, INTENT(IN)                                :: conjugator, max_iter
     731             :       REAL(KIND=dp), INTENT(IN)                          :: eps_convergence, eps_filter
     732             :       LOGICAL, INTENT(OUT)                               :: converged
     733             : 
     734             :       CHARACTER(len=*), PARAMETER :: routineN = 'solve_riccati_equation'
     735             : 
     736             :       INTEGER                                            :: handle, istep, iteration, nsteps, &
     737             :                                                             unit_nr, update_prec_freq
     738             :       LOGICAL                                            :: prepare_to_exit, present_oo, present_vv, &
     739             :                                                             quadratic_term, restart_conjugator
     740             :       REAL(KIND=dp)                                      :: best_norm, best_step_size, beta, c0, c1, &
     741             :                                                             c2, c3, denom, kappa, numer, &
     742             :                                                             obj_function, t1, t2, tau
     743             :       REAL(KIND=dp), DIMENSION(3)                        :: step_size
     744             :       TYPE(cp_logger_type), POINTER                      :: logger
     745             :       TYPE(dbcsr_type)                                   :: aux1, aux2, grad, m, n, oo1, oo2, prec, &
     746             :                                                             res_trial, step, step_oo, vv_step
     747             : 
     748             : !TYPE(dbcsr_type)                      :: qqqq, pppp, zero_pq, zero_qp
     749             : 
     750           0 :       CALL timeset(routineN, handle)
     751             : 
     752           0 :       logger => cp_get_default_logger()
     753           0 :       IF (logger%para_env%is_source()) THEN
     754           0 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     755             :       ELSE
     756             :          unit_nr = -1
     757             :       END IF
     758             : 
     759           0 :       t1 = m_walltime()
     760             : 
     761             : !IF (level.gt.5) THEN
     762             : !  CPErrorMessage(cp_failure_level,routineP,"recursion level is too high")
     763             : !  CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
     764             : !ENDIF
     765             : !IF (unit_nr>0) THEN
     766             : !   WRITE(unit_nr,*) &
     767             : !      "========== LEVEL ",level,"=========="
     768             : !ENDIF
     769             : !CALL dbcsr_print(qq)
     770             : !CALL dbcsr_print(pp)
     771             : !CALL dbcsr_print(qp)
     772             : !!CALL dbcsr_print(pq)
     773             : !IF (unit_nr>0) THEN
     774             : !   WRITE(unit_nr,*) &
     775             : !      "====== END LEVEL ",level,"=========="
     776             : !ENDIF
     777             : 
     778           0 :       quadratic_term = .NOT. neglect_quadratic_term
     779           0 :       present_oo = PRESENT(oo)
     780           0 :       present_vv = PRESENT(vv)
     781             : 
     782             :       ! create aux1 matrix and init
     783           0 :       CALL dbcsr_create(aux1, template=pp)
     784           0 :       CALL dbcsr_copy(aux1, pp)
     785           0 :       CALL dbcsr_scale(aux1, -1.0_dp)
     786             : 
     787             :       ! create aux2 matrix and init
     788           0 :       CALL dbcsr_create(aux2, template=qq)
     789           0 :       CALL dbcsr_copy(aux2, qq)
     790             : 
     791             :       ! create the gradient matrix and init
     792           0 :       CALL dbcsr_create(grad, template=x)
     793           0 :       CALL dbcsr_set(grad, 0.0_dp)
     794             : 
     795             :       ! create a preconditioner
     796             :       ! RZK-warning how to apply it to up_down tensor?
     797           0 :       CALL dbcsr_create(prec, template=x)
     798             :       !CALL create_preconditioner(prec,aux1,aux2,qp,res,tensor_type,eps_filter)
     799             :       !CALL dbcsr_set(prec,1.0_dp)
     800             : 
     801             :       ! create the step matrix and init
     802           0 :       CALL dbcsr_create(step, template=x)
     803             :       !CALL dbcsr_hadamard_product(prec,grad,step)
     804             :       !CALL dbcsr_scale(step,-1.0_dp)
     805             : 
     806           0 :       CALL dbcsr_create(n, template=x)
     807           0 :       CALL dbcsr_create(m, template=x)
     808           0 :       CALL dbcsr_create(oo1, template=pp)
     809           0 :       CALL dbcsr_create(oo2, template=pp)
     810           0 :       CALL dbcsr_create(res_trial, template=res)
     811           0 :       CALL dbcsr_create(vv_step, template=res)
     812           0 :       CALL dbcsr_create(step_oo, template=res)
     813             : 
     814             :       ! start conjugate gradient iterations
     815           0 :       iteration = 0
     816           0 :       converged = .FALSE.
     817           0 :       prepare_to_exit = .FALSE.
     818           0 :       beta = 0.0_dp
     819           0 :       best_step_size = 0.0_dp
     820           0 :       best_norm = 1.0E+100_dp
     821             :       !ecorr=0.0_dp
     822             :       !change_ecorr=0.0_dp
     823           0 :       restart_conjugator = .FALSE.
     824           0 :       update_prec_freq = 20
     825             :       DO
     826             : 
     827             :          ! (re)-compute the residuals
     828           0 :          IF (iteration .EQ. 0) THEN
     829           0 :             CALL dbcsr_copy(res, qp)
     830           0 :             IF (present_oo) THEN
     831             :                CALL dbcsr_multiply("N", "N", +1.0_dp, qq, x, 0.0_dp, res_trial, &
     832           0 :                                    filter_eps=eps_filter)
     833             :                CALL dbcsr_multiply("N", "N", +1.0_dp, res_trial, oo, 1.0_dp, res, &
     834           0 :                                    filter_eps=eps_filter)
     835             :             ELSE
     836             :                CALL dbcsr_multiply("N", "N", +1.0_dp, qq, x, 1.0_dp, res, &
     837           0 :                                    filter_eps=eps_filter)
     838             :             END IF
     839           0 :             IF (present_vv) THEN
     840             :                CALL dbcsr_multiply("N", "N", -1.0_dp, x, pp, 0.0_dp, res_trial, &
     841           0 :                                    filter_eps=eps_filter)
     842             :                CALL dbcsr_multiply("N", "N", +1.0_dp, vv, res_trial, 1.0_dp, res, &
     843           0 :                                    filter_eps=eps_filter)
     844             :             ELSE
     845             :                CALL dbcsr_multiply("N", "N", -1.0_dp, x, pp, 1.0_dp, res, &
     846           0 :                                    filter_eps=eps_filter)
     847             :             END IF
     848           0 :             IF (quadratic_term) THEN
     849           0 :                IF (present_oo) THEN
     850             :                   CALL dbcsr_multiply("N", "N", +1.0_dp, pq, x, 0.0_dp, oo1, &
     851           0 :                                       filter_eps=eps_filter)
     852             :                   CALL dbcsr_multiply("N", "N", +1.0_dp, oo1, oo, 0.0_dp, oo2, &
     853           0 :                                       filter_eps=eps_filter)
     854             :                ELSE
     855             :                   CALL dbcsr_multiply("N", "N", +1.0_dp, pq, x, 0.0_dp, oo2, &
     856           0 :                                       filter_eps=eps_filter)
     857             :                END IF
     858           0 :                IF (present_vv) THEN
     859             :                   CALL dbcsr_multiply("N", "N", -1.0_dp, x, oo2, 0.0_dp, res_trial, &
     860           0 :                                       filter_eps=eps_filter)
     861             :                   CALL dbcsr_multiply("N", "N", +1.0_dp, vv, res_trial, 1.0_dp, res, &
     862           0 :                                       filter_eps=eps_filter)
     863             :                ELSE
     864             :                   CALL dbcsr_multiply("N", "N", -1.0_dp, x, oo2, 1.0_dp, res, &
     865           0 :                                       filter_eps=eps_filter)
     866             :                END IF
     867             :             END IF
     868           0 :             best_norm = dbcsr_maxabs(res)
     869             :          ELSE
     870           0 :             CALL dbcsr_add(res, m, 1.0_dp, best_step_size)
     871           0 :             CALL dbcsr_add(res, n, 1.0_dp, -best_step_size*best_step_size)
     872           0 :             CALL dbcsr_filter(res, eps_filter)
     873             :          END IF
     874             : 
     875             :          ! check convergence and other exit criteria
     876           0 :          converged = (best_norm .LT. eps_convergence)
     877           0 :          IF (converged .OR. (iteration .GE. max_iter)) THEN
     878             :             prepare_to_exit = .TRUE.
     879             :          END IF
     880             : 
     881           0 :          IF (.NOT. prepare_to_exit) THEN
     882             : 
     883             :             ! update aux1=-pp-pq.x.oo and aux2=qq-vv.x.pq
     884           0 :             IF (quadratic_term) THEN
     885           0 :                IF (iteration .EQ. 0) THEN
     886           0 :                   IF (present_oo) THEN
     887             :                      CALL dbcsr_multiply("N", "N", -1.0_dp, pq, x, 0.0_dp, oo1, &
     888           0 :                                          filter_eps=eps_filter)
     889             :                      CALL dbcsr_multiply("N", "N", +1.0_dp, oo1, oo, 1.0_dp, aux1, &
     890           0 :                                          filter_eps=eps_filter)
     891             :                   ELSE
     892             :                      CALL dbcsr_multiply("N", "N", -1.0_dp, pq, x, 1.0_dp, aux1, &
     893           0 :                                          filter_eps=eps_filter)
     894             :                   END IF
     895           0 :                   IF (present_vv) THEN
     896             :                      CALL dbcsr_multiply("N", "N", -1.0_dp, vv, x, 0.0_dp, res_trial, &
     897           0 :                                          filter_eps=eps_filter)
     898             :                      CALL dbcsr_multiply("N", "N", +1.0_dp, res_trial, pq, 1.0_dp, aux2, &
     899           0 :                                          filter_eps=eps_filter)
     900             :                   ELSE
     901             :                      CALL dbcsr_multiply("N", "N", -1.0_dp, x, pq, 1.0_dp, aux2, &
     902           0 :                                          filter_eps=eps_filter)
     903             :                   END IF
     904             :                ELSE
     905           0 :                   IF (present_oo) THEN
     906             :                      CALL dbcsr_multiply("N", "N", -best_step_size, pq, step_oo, 1.0_dp, aux1, &
     907           0 :                                          filter_eps=eps_filter)
     908             :                   ELSE
     909             :                      CALL dbcsr_multiply("N", "N", -best_step_size, pq, step, 1.0_dp, aux1, &
     910           0 :                                          filter_eps=eps_filter)
     911             :                   END IF
     912           0 :                   IF (present_vv) THEN
     913             :                      CALL dbcsr_multiply("N", "N", -best_step_size, vv_step, pq, 1.0_dp, aux2, &
     914           0 :                                          filter_eps=eps_filter)
     915             :                   ELSE
     916             :                      CALL dbcsr_multiply("N", "N", -best_step_size, step, pq, 1.0_dp, aux2, &
     917           0 :                                          filter_eps=eps_filter)
     918             :                   END IF
     919             :                END IF
     920             :             END IF
     921             : 
     922             :             ! recompute the gradient, do not update it yet
     923             :             ! use m matrix as a temporary storage
     924             :             ! grad=t(vv).res.t(aux1)+t(aux2).res.t(oo)
     925           0 :             IF (present_vv) THEN
     926             :                CALL dbcsr_multiply("N", "T", 1.0_dp, res, aux1, 0.0_dp, res_trial, &
     927           0 :                                    filter_eps=eps_filter)
     928             :                CALL dbcsr_multiply("T", "N", 1.0_dp, vv, res_trial, 0.0_dp, m, &
     929           0 :                                    filter_eps=eps_filter)
     930             :             ELSE
     931             :                CALL dbcsr_multiply("N", "T", 1.0_dp, res, aux1, 0.0_dp, m, &
     932           0 :                                    filter_eps=eps_filter)
     933             :             END IF
     934           0 :             IF (present_oo) THEN
     935             :                CALL dbcsr_multiply("T", "N", 1.0_dp, aux1, res, 0.0_dp, res_trial, &
     936           0 :                                    filter_eps=eps_filter)
     937             :                CALL dbcsr_multiply("N", "T", 1.0_dp, res_trial, oo, 1.0_dp, m, &
     938           0 :                                    filter_eps=eps_filter)
     939             :             ELSE
     940             :                CALL dbcsr_multiply("T", "N", 1.0_dp, aux2, res, 1.0_dp, m, &
     941           0 :                                    filter_eps=eps_filter)
     942             :             END IF
     943             : 
     944             :             ! compute preconditioner
     945             :             !IF (iteration.eq.0.OR.(mod(iteration,update_prec_freq).eq.0)) THEN
     946           0 :             IF (iteration .EQ. 0) THEN
     947           0 :                CALL create_preconditioner(prec, aux1, aux2, eps_filter)
     948             :                !restart_conjugator=.TRUE.
     949             : !CALL dbcsr_set(prec,1.0_dp)
     950             : !CALL dbcsr_print(prec)
     951             :             END IF
     952             : 
     953             :             ! compute the conjugation coefficient - beta
     954           0 :             IF ((iteration .EQ. 0) .OR. restart_conjugator) THEN
     955           0 :                beta = 0.0_dp
     956             :             ELSE
     957           0 :                restart_conjugator = .FALSE.
     958           0 :                SELECT CASE (conjugator)
     959             :                CASE (cg_hestenes_stiefel)
     960           0 :                   CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
     961           0 :                   CALL dbcsr_hadamard_product(prec, grad, n)
     962           0 :                   CALL dbcsr_dot(n, m, numer)
     963           0 :                   CALL dbcsr_dot(grad, step, denom)
     964           0 :                   beta = numer/denom
     965             :                CASE (cg_fletcher_reeves)
     966           0 :                   CALL dbcsr_hadamard_product(prec, grad, n)
     967           0 :                   CALL dbcsr_dot(grad, n, denom)
     968           0 :                   CALL dbcsr_hadamard_product(prec, m, n)
     969           0 :                   CALL dbcsr_dot(m, n, numer)
     970           0 :                   beta = numer/denom
     971             :                CASE (cg_polak_ribiere)
     972           0 :                   CALL dbcsr_hadamard_product(prec, grad, n)
     973           0 :                   CALL dbcsr_dot(grad, n, denom)
     974           0 :                   CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
     975           0 :                   CALL dbcsr_hadamard_product(prec, grad, n)
     976           0 :                   CALL dbcsr_dot(n, m, numer)
     977           0 :                   beta = numer/denom
     978             :                CASE (cg_fletcher)
     979           0 :                   CALL dbcsr_hadamard_product(prec, m, n)
     980           0 :                   CALL dbcsr_dot(m, n, numer)
     981           0 :                   CALL dbcsr_dot(grad, step, denom)
     982           0 :                   beta = -1.0_dp*numer/denom
     983             :                CASE (cg_liu_storey)
     984           0 :                   CALL dbcsr_dot(grad, step, denom)
     985           0 :                   CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
     986           0 :                   CALL dbcsr_hadamard_product(prec, grad, n)
     987           0 :                   CALL dbcsr_dot(n, m, numer)
     988           0 :                   beta = -1.0_dp*numer/denom
     989             :                CASE (cg_dai_yuan)
     990           0 :                   CALL dbcsr_hadamard_product(prec, m, n)
     991           0 :                   CALL dbcsr_dot(m, n, numer)
     992           0 :                   CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
     993           0 :                   CALL dbcsr_dot(grad, step, denom)
     994           0 :                   beta = numer/denom
     995             :                CASE (cg_hager_zhang)
     996           0 :                   CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
     997           0 :                   CALL dbcsr_dot(grad, step, denom)
     998           0 :                   CALL dbcsr_hadamard_product(prec, grad, n)
     999           0 :                   CALL dbcsr_dot(n, grad, numer)
    1000           0 :                   kappa = 2.0_dp*numer/denom
    1001           0 :                   CALL dbcsr_dot(n, m, numer)
    1002           0 :                   tau = numer/denom
    1003           0 :                   CALL dbcsr_dot(step, m, numer)
    1004           0 :                   beta = tau - kappa*numer/denom
    1005             :                CASE (cg_zero)
    1006           0 :                   beta = 0.0_dp
    1007             :                CASE DEFAULT
    1008           0 :                   CPABORT("illegal conjugator")
    1009             :                END SELECT
    1010             :             END IF ! iteration.eq.0
    1011             : 
    1012             :             ! move the current gradient to its storage
    1013           0 :             CALL dbcsr_copy(grad, m)
    1014             : 
    1015             :             ! precondition new gradient (use m as tmp storage)
    1016           0 :             CALL dbcsr_hadamard_product(prec, grad, m)
    1017           0 :             CALL dbcsr_filter(m, eps_filter)
    1018             : 
    1019             :             ! recompute the step direction
    1020           0 :             CALL dbcsr_add(step, m, beta, -1.0_dp)
    1021           0 :             CALL dbcsr_filter(step, eps_filter)
    1022             : 
    1023             : !! ALTERNATIVE METHOD TO OBTAIN THE STEP FROM THE GRADIENT
    1024             : !CALL dbcsr_init(qqqq)
    1025             : !CALL dbcsr_create(qqqq,template=qq)
    1026             : !CALL dbcsr_init(pppp)
    1027             : !CALL dbcsr_create(pppp,template=pp)
    1028             : !CALL dbcsr_init(zero_pq)
    1029             : !CALL dbcsr_create(zero_pq,template=pq)
    1030             : !CALL dbcsr_init(zero_qp)
    1031             : !CALL dbcsr_create(zero_qp,template=qp)
    1032             : !CALL dbcsr_multiply("T","N",1.0_dp,aux2,aux2,0.0_dp,qqqq,&
    1033             : !        filter_eps=eps_filter)
    1034             : !CALL dbcsr_multiply("N","T",-1.0_dp,aux1,aux1,0.0_dp,pppp,&
    1035             : !        filter_eps=eps_filter)
    1036             : !CALL dbcsr_set(zero_qp,0.0_dp)
    1037             : !CALL dbcsr_set(zero_pq,0.0_dp)
    1038             : !CALL solve_riccati_equation(pppp,qqqq,grad,zero_pq,zero_qp,zero_qp,&
    1039             : !               .TRUE.,tensor_type,&
    1040             : !               conjugator,max_iter,eps_convergence,eps_filter,&
    1041             : !               converged,level+1)
    1042             : !CALL dbcsr_release(qqqq)
    1043             : !CALL dbcsr_release(pppp)
    1044             : !CALL dbcsr_release(zero_qp)
    1045             : !CALL dbcsr_release(zero_pq)
    1046             : 
    1047             :             ! calculate the optimal step size
    1048             :             ! m=step.aux1+aux2.step
    1049           0 :             IF (present_vv) THEN
    1050             :                CALL dbcsr_multiply("N", "N", 1.0_dp, vv, step, 0.0_dp, vv_step, &
    1051           0 :                                    filter_eps=eps_filter)
    1052             :                CALL dbcsr_multiply("N", "N", 1.0_dp, vv_step, aux1, 0.0_dp, m, &
    1053           0 :                                    filter_eps=eps_filter)
    1054             :             ELSE
    1055             :                CALL dbcsr_multiply("N", "N", 1.0_dp, step, aux1, 0.0_dp, m, &
    1056           0 :                                    filter_eps=eps_filter)
    1057             :             END IF
    1058           0 :             IF (present_oo) THEN
    1059             :                CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo, 0.0_dp, step_oo, &
    1060           0 :                                    filter_eps=eps_filter)
    1061             :                CALL dbcsr_multiply("N", "N", 1.0_dp, aux2, step_oo, 1.0_dp, m, &
    1062           0 :                                    filter_eps=eps_filter)
    1063             :             ELSE
    1064             :                CALL dbcsr_multiply("N", "N", 1.0_dp, aux2, step, 1.0_dp, m, &
    1065           0 :                                    filter_eps=eps_filter)
    1066             :             END IF
    1067             : 
    1068           0 :             IF (quadratic_term) THEN
    1069             :                ! n=step.pq.step
    1070           0 :                IF (present_oo) THEN
    1071             :                   CALL dbcsr_multiply("N", "N", 1.0_dp, pq, step, 0.0_dp, oo1, &
    1072           0 :                                       filter_eps=eps_filter)
    1073             :                   CALL dbcsr_multiply("N", "N", 1.0_dp, oo1, oo, 0.0_dp, oo2, &
    1074           0 :                                       filter_eps=eps_filter)
    1075             :                ELSE
    1076             :                   CALL dbcsr_multiply("N", "N", 1.0_dp, pq, step, 0.0_dp, oo2, &
    1077           0 :                                       filter_eps=eps_filter)
    1078             :                END IF
    1079           0 :                IF (present_vv) THEN
    1080             :                   CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo2, 0.0_dp, res_trial, &
    1081           0 :                                       filter_eps=eps_filter)
    1082             :                   CALL dbcsr_multiply("N", "N", 1.0_dp, vv, res_trial, 0.0_dp, n, &
    1083           0 :                                       filter_eps=eps_filter)
    1084             :                ELSE
    1085             :                   CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo2, 0.0_dp, n, &
    1086           0 :                                       filter_eps=eps_filter)
    1087             :                END IF
    1088             : 
    1089             :             ELSE
    1090           0 :                CALL dbcsr_set(n, 0.0_dp)
    1091             :             END IF
    1092             : 
    1093             :             ! calculate coefficients of the cubic eq for alpha - step size
    1094           0 :             c0 = 2.0_dp*(dbcsr_frobenius_norm(n))**2
    1095             : 
    1096           0 :             CALL dbcsr_dot(m, n, c1)
    1097           0 :             c1 = -3.0_dp*c1
    1098             : 
    1099           0 :             CALL dbcsr_dot(res, n, c2)
    1100           0 :             c2 = -2.0_dp*c2 + (dbcsr_frobenius_norm(m))**2
    1101             : 
    1102           0 :             CALL dbcsr_dot(res, m, c3)
    1103             : 
    1104             :             ! find step size
    1105           0 :             CALL analytic_line_search(c0, c1, c2, c3, step_size, nsteps)
    1106             : 
    1107           0 :             IF (nsteps .EQ. 0) THEN
    1108           0 :                CPABORT("no step sizes!")
    1109             :             END IF
    1110             :             ! if we have several possible step sizes
    1111             :             ! choose one with the lowest objective function
    1112           0 :             best_norm = 1.0E+100_dp
    1113           0 :             best_step_size = 0.0_dp
    1114           0 :             DO istep = 1, nsteps
    1115             :                ! recompute the residues
    1116           0 :                CALL dbcsr_copy(res_trial, res)
    1117           0 :                CALL dbcsr_add(res_trial, m, 1.0_dp, step_size(istep))
    1118           0 :                CALL dbcsr_add(res_trial, n, 1.0_dp, -step_size(istep)*step_size(istep))
    1119           0 :                CALL dbcsr_filter(res_trial, eps_filter)
    1120             :                ! RZK-warning objective function might be different in the case of
    1121             :                ! tensor_up_down
    1122             :                !obj_function=0.5_dp*(dbcsr_frobenius_norm(res_trial))**2
    1123           0 :                obj_function = dbcsr_maxabs(res_trial)
    1124           0 :                IF (obj_function .LT. best_norm) THEN
    1125           0 :                   best_norm = obj_function
    1126           0 :                   best_step_size = step_size(istep)
    1127             :                END IF
    1128             :             END DO
    1129             : 
    1130             :          END IF
    1131             : 
    1132             :          ! update X along the line
    1133           0 :          CALL dbcsr_add(x, step, 1.0_dp, best_step_size)
    1134           0 :          CALL dbcsr_filter(x, eps_filter)
    1135             : 
    1136             :          ! evaluate current energy correction
    1137             :          !change_ecorr=ecorr
    1138             :          !CALL dbcsr_dot(qp,x,ecorr,"T","N")
    1139             :          !change_ecorr=ecorr-change_ecorr
    1140             : 
    1141             :          ! check convergence and other exit criteria
    1142           0 :          converged = (best_norm .LT. eps_convergence)
    1143           0 :          IF (converged .OR. (iteration .GE. max_iter)) THEN
    1144           0 :             prepare_to_exit = .TRUE.
    1145             :          END IF
    1146             : 
    1147           0 :          t2 = m_walltime()
    1148             : 
    1149           0 :          IF (unit_nr > 0) THEN
    1150             :             WRITE (unit_nr, '(T6,A,1X,I4,1X,E12.3,F8.3)') &
    1151           0 :                "RICCATI iter ", iteration, best_norm, t2 - t1
    1152             :             !WRITE(unit_nr,'(T6,A,1X,I4,1X,F15.9,F15.9,E12.3,F8.3)') &
    1153             :             !   "RICCATI iter ",iteration,ecorr,change_ecorr,best_norm,t2-t1
    1154             :          END IF
    1155             : 
    1156           0 :          t1 = m_walltime()
    1157             : 
    1158           0 :          iteration = iteration + 1
    1159             : 
    1160           0 :          IF (prepare_to_exit) EXIT
    1161             : 
    1162             :       END DO
    1163             : 
    1164           0 :       CALL dbcsr_release(aux1)
    1165           0 :       CALL dbcsr_release(aux2)
    1166           0 :       CALL dbcsr_release(grad)
    1167           0 :       CALL dbcsr_release(step)
    1168           0 :       CALL dbcsr_release(n)
    1169           0 :       CALL dbcsr_release(m)
    1170           0 :       CALL dbcsr_release(oo1)
    1171           0 :       CALL dbcsr_release(oo2)
    1172           0 :       CALL dbcsr_release(res_trial)
    1173           0 :       CALL dbcsr_release(vv_step)
    1174           0 :       CALL dbcsr_release(step_oo)
    1175             : 
    1176           0 :       CALL timestop(handle)
    1177             : 
    1178           0 :    END SUBROUTINE solve_riccati_equation
    1179             : 
    1180             : ! **************************************************************************************************
    1181             : !> \brief Computes a preconditioner from diagonal elements of ~f_oo, ~f_vv
    1182             : !>        The preconditioner is approximately equal to
    1183             : !>        prec_ai ~ (e_a - e_i)^(-2)
    1184             : !>        However, the real expression is more complex
    1185             : !> \param prec ...
    1186             : !> \param pp ...
    1187             : !> \param qq ...
    1188             : !> \param eps_filter ...
    1189             : !> \par History
    1190             : !>       2011.07 created [Rustam Z Khaliullin]
    1191             : !> \author Rustam Z Khaliullin
    1192             : ! **************************************************************************************************
    1193           0 :    SUBROUTINE create_preconditioner(prec, pp, qq, eps_filter)
    1194             : 
    1195             :       TYPE(dbcsr_type), INTENT(OUT)                      :: prec
    1196             :       TYPE(dbcsr_type), INTENT(IN)                       :: pp, qq
    1197             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
    1198             : 
    1199             :       CHARACTER(len=*), PARAMETER :: routineN = 'create_preconditioner'
    1200             : 
    1201             :       INTEGER                                            :: col, handle, hold, iblock_col, &
    1202             :                                                             iblock_row, mynode, nblkcols_tot, &
    1203             :                                                             nblkrows_tot, p_nrows, q_nrows, row
    1204             :       LOGICAL                                            :: tr
    1205           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: p_diagonal, q_diagonal
    1206           0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_new_block
    1207             :       TYPE(dbcsr_distribution_type)                      :: dist
    1208             :       TYPE(dbcsr_type)                                   :: pp_diag, qq_diag, t1, t2, tmp
    1209             : 
    1210             : !LOGICAL, INTENT(IN)                      :: use_virt_orbs
    1211             : 
    1212           0 :       CALL timeset(routineN, handle)
    1213             : 
    1214             : !    ! copy diagonal elements
    1215             : !    CALL dbcsr_get_info(pp,nfullrows_total=nrows)
    1216             : !    CALL dbcsr_init(pp_diag)
    1217             : !    CALL dbcsr_create(pp_diag,template=pp)
    1218             : !    ALLOCATE(diagonal(nrows))
    1219             : !    CALL dbcsr_get_diag(pp,diagonal)
    1220             : !    CALL dbcsr_add_on_diag(pp_diag,1.0_dp)
    1221             : !    CALL dbcsr_set_diag(pp_diag,diagonal)
    1222             : !    DEALLOCATE(diagonal)
    1223             : !
    1224             :       ! initialize a matrix to 1.0
    1225           0 :       CALL dbcsr_create(tmp, template=prec)
    1226             :       ! use an ugly hack to set all elements of tmp to 1
    1227             :       ! because dbcsr_set does not do it (despite its name)
    1228             :       !CALL dbcsr_set(tmp,1.0_dp)
    1229           0 :       CALL dbcsr_get_info(tmp, distribution=dist)
    1230           0 :       CALL dbcsr_distribution_get(dist, mynode=mynode)
    1231           0 :       CALL dbcsr_work_create(tmp, work_mutable=.TRUE.)
    1232           0 :       nblkrows_tot = dbcsr_nblkrows_total(tmp)
    1233           0 :       nblkcols_tot = dbcsr_nblkcols_total(tmp)
    1234           0 :       DO row = 1, nblkrows_tot
    1235           0 :          DO col = 1, nblkcols_tot
    1236           0 :             tr = .FALSE.
    1237           0 :             iblock_row = row
    1238           0 :             iblock_col = col
    1239           0 :             CALL dbcsr_get_stored_coordinates(tmp, iblock_row, iblock_col, hold)
    1240           0 :             IF (hold .EQ. mynode) THEN
    1241           0 :                NULLIFY (p_new_block)
    1242           0 :                CALL dbcsr_reserve_block2d(tmp, iblock_row, iblock_col, p_new_block)
    1243           0 :                CPASSERT(ASSOCIATED(p_new_block))
    1244           0 :                p_new_block(:, :) = 1.0_dp
    1245             :             END IF ! mynode
    1246             :          END DO
    1247             :       END DO
    1248           0 :       CALL dbcsr_finalize(tmp)
    1249             : 
    1250             :       ! copy diagonal elements of pp into cols of a matrix
    1251           0 :       CALL dbcsr_get_info(pp, nfullrows_total=p_nrows)
    1252           0 :       CALL dbcsr_create(pp_diag, template=pp)
    1253           0 :       ALLOCATE (p_diagonal(p_nrows))
    1254           0 :       CALL dbcsr_get_diag(pp, p_diagonal)
    1255           0 :       CALL dbcsr_add_on_diag(pp_diag, 1.0_dp)
    1256           0 :       CALL dbcsr_set_diag(pp_diag, p_diagonal)
    1257             :       ! RZK-warning is it possible to use dbcsr_scale_by_vector?
    1258             :       ! or even insert elements directly in the prev cycles
    1259           0 :       CALL dbcsr_create(t2, template=prec)
    1260             :       CALL dbcsr_multiply("N", "N", 1.0_dp, tmp, pp_diag, &
    1261           0 :                           0.0_dp, t2, filter_eps=eps_filter)
    1262             : 
    1263             :       ! copy diagonal elements qq into rows of a matrix
    1264           0 :       CALL dbcsr_get_info(qq, nfullrows_total=q_nrows)
    1265           0 :       CALL dbcsr_create(qq_diag, template=qq)
    1266           0 :       ALLOCATE (q_diagonal(q_nrows))
    1267           0 :       CALL dbcsr_get_diag(qq, q_diagonal)
    1268           0 :       CALL dbcsr_add_on_diag(qq_diag, 1.0_dp)
    1269           0 :       CALL dbcsr_set_diag(qq_diag, q_diagonal)
    1270           0 :       CALL dbcsr_set(tmp, 1.0_dp)
    1271           0 :       CALL dbcsr_create(t1, template=prec)
    1272             :       CALL dbcsr_multiply("N", "N", 1.0_dp, qq_diag, tmp, &
    1273           0 :                           0.0_dp, t1, filter_eps=eps_filter)
    1274             : 
    1275           0 :       CALL dbcsr_hadamard_product(t1, t2, prec)
    1276           0 :       CALL dbcsr_release(t1)
    1277           0 :       CALL dbcsr_scale(prec, 2.0_dp)
    1278             : 
    1279             :       ! Get the diagonal of tr(qq).qq
    1280             :       CALL dbcsr_multiply("T", "N", 1.0_dp, qq, qq, &
    1281             :                           0.0_dp, qq_diag, retain_sparsity=.TRUE., &
    1282           0 :                           filter_eps=eps_filter)
    1283           0 :       CALL dbcsr_get_diag(qq_diag, q_diagonal)
    1284           0 :       CALL dbcsr_set(qq_diag, 0.0_dp)
    1285           0 :       CALL dbcsr_add_on_diag(qq_diag, 1.0_dp)
    1286           0 :       CALL dbcsr_set_diag(qq_diag, q_diagonal)
    1287           0 :       DEALLOCATE (q_diagonal)
    1288           0 :       CALL dbcsr_set(tmp, 1.0_dp)
    1289             :       CALL dbcsr_multiply("N", "N", 1.0_dp, qq_diag, tmp, &
    1290           0 :                           0.0_dp, t2, filter_eps=eps_filter)
    1291           0 :       CALL dbcsr_release(qq_diag)
    1292           0 :       CALL dbcsr_add(prec, t2, 1.0_dp, 1.0_dp)
    1293             : 
    1294             :       ! Get the diagonal of pp.tr(pp)
    1295             :       CALL dbcsr_multiply("N", "T", 1.0_dp, pp, pp, &
    1296             :                           0.0_dp, pp_diag, retain_sparsity=.TRUE., &
    1297           0 :                           filter_eps=eps_filter)
    1298           0 :       CALL dbcsr_get_diag(pp_diag, p_diagonal)
    1299           0 :       CALL dbcsr_set(pp_diag, 0.0_dp)
    1300           0 :       CALL dbcsr_add_on_diag(pp_diag, 1.0_dp)
    1301           0 :       CALL dbcsr_set_diag(pp_diag, p_diagonal)
    1302           0 :       DEALLOCATE (p_diagonal)
    1303           0 :       CALL dbcsr_set(tmp, 1.0_dp)
    1304             :       CALL dbcsr_multiply("N", "N", 1.0_dp, tmp, pp_diag, &
    1305           0 :                           0.0_dp, t2, filter_eps=eps_filter)
    1306           0 :       CALL dbcsr_release(tmp)
    1307           0 :       CALL dbcsr_release(pp_diag)
    1308           0 :       CALL dbcsr_add(prec, t2, 1.0_dp, 1.0_dp)
    1309             : 
    1310             :       ! now add the residual component
    1311             :       !CALL dbcsr_hadamard_product(res,qp,t2)
    1312             :       !CALL dbcsr_add(prec,t2,1.0_dp,-2.0_dp)
    1313           0 :       CALL dbcsr_release(t2)
    1314           0 :       CALL inverse_of_elements(prec)
    1315           0 :       CALL dbcsr_filter(prec, eps_filter)
    1316             : 
    1317           0 :       CALL timestop(handle)
    1318             : 
    1319           0 :    END SUBROUTINE create_preconditioner
    1320             : 
    1321             : ! **************************************************************************************************
    1322             : !> \brief Computes 1/x of the matrix elements.
    1323             : !> \param matrix ...
    1324             : !> \author Ole Schuett
    1325             : ! **************************************************************************************************
    1326           0 :    SUBROUTINE inverse_of_elements(matrix)
    1327             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
    1328             : 
    1329             :       CHARACTER(len=*), PARAMETER :: routineN = 'inverse_of_elements'
    1330             : 
    1331             :       INTEGER                                            :: handle
    1332           0 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: block
    1333             :       TYPE(dbcsr_iterator_type)                          :: iter
    1334             : 
    1335           0 :       CALL timeset(routineN, handle)
    1336           0 :       CALL dbcsr_iterator_start(iter, matrix)
    1337           0 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
    1338           0 :          CALL dbcsr_iterator_next_block(iter, block=block)
    1339           0 :          block = 1.0_dp/block
    1340             :       END DO
    1341           0 :       CALL dbcsr_iterator_stop(iter)
    1342           0 :       CALL timestop(handle)
    1343             : 
    1344           0 :    END SUBROUTINE inverse_of_elements
    1345             : 
    1346             : ! **************************************************************************************************
    1347             : !> \brief Finds real roots of a cubic equation
    1348             : !>    >        a*x**3 + b*x**2 + c*x + d = 0
    1349             : !>        and returns only those roots for which the derivative is positive
    1350             : !>
    1351             : !>   Step 0: Check the true order of the equation. Cubic, quadratic, linear?
    1352             : !>   Step 1: Calculate p and q
    1353             : !>           p = ( 3*c/a - (b/a)**2 ) / 3
    1354             : !>           q = ( 2*(b/a)**3 - 9*b*c/a/a + 27*d/a ) / 27
    1355             : !>   Step 2: Calculate discriminant D
    1356             : !>           D = (p/3)**3 + (q/2)**2
    1357             : !>   Step 3: Depending on the sign of D, we follow different strategy.
    1358             : !>           If D<0, three distinct real roots.
    1359             : !>           If D=0, three real roots of which at least two are equal.
    1360             : !>           If D>0, one real and two complex roots.
    1361             : !>   Step 3a: For D>0 and D=0,
    1362             : !>           Calculate u and v
    1363             : !>           u = cubic_root(-q/2 + sqrt(D))
    1364             : !>           v = cubic_root(-q/2 - sqrt(D))
    1365             : !>           Find the three transformed roots
    1366             : !>           y1 = u + v
    1367             : !>           y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2
    1368             : !>           y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2
    1369             : !>   Step 3b Alternately, for D<0, a trigonometric formulation is more convenient
    1370             : !>           y1 =  2 * sqrt(|p|/3) * cos(phi/3)
    1371             : !>           y2 = -2 * sqrt(|p|/3) * cos((phi+pi)/3)
    1372             : !>           y3 = -2 * sqrt(|p|/3) * cos((phi-pi)/3)
    1373             : !>           where phi = acos(-q/2/sqrt(|p|**3/27))
    1374             : !>                 pi  = 3.141592654...
    1375             : !>   Step 4  Find the real roots
    1376             : !>           x = y - b/a/3
    1377             : !>   Step 5  Check the derivative and return only those real roots
    1378             : !>           for which the derivative is positive
    1379             : !>
    1380             : !> \param a ...
    1381             : !> \param b ...
    1382             : !> \param c ...
    1383             : !> \param d ...
    1384             : !> \param minima ...
    1385             : !> \param nmins ...
    1386             : !> \par History
    1387             : !>       2011.06 created [Rustam Z Khaliullin]
    1388             : !> \author Rustam Z Khaliullin
    1389             : ! **************************************************************************************************
    1390           0 :    SUBROUTINE analytic_line_search(a, b, c, d, minima, nmins)
    1391             : 
    1392             :       REAL(KIND=dp), INTENT(IN)                          :: a, b, c, d
    1393             :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT)           :: minima
    1394             :       INTEGER, INTENT(OUT)                               :: nmins
    1395             : 
    1396             :       INTEGER                                            :: i, nroots
    1397             :       REAL(KIND=dp)                                      :: DD, der, p, phi, pi, q, temp1, temp2, u, &
    1398             :                                                             v, y1, y2, y2i, y2r, y3
    1399             :       REAL(KIND=dp), DIMENSION(3)                        :: x
    1400             : 
    1401             : !    CALL timeset(routineN,handle)
    1402             : 
    1403           0 :       pi = ACOS(-1.0_dp)
    1404             : 
    1405             :       ! Step 0: Check coefficients and find the true order of the eq
    1406           0 :       IF (a .EQ. 0.0_dp) THEN
    1407           0 :          IF (b .EQ. 0.0_dp) THEN
    1408           0 :             IF (c .EQ. 0.0_dp) THEN
    1409             :                ! Non-equation, no valid solutions
    1410             :                nroots = 0
    1411             :             ELSE
    1412             :                ! Linear equation with one root.
    1413           0 :                nroots = 1
    1414           0 :                x(1) = -d/c
    1415             :             END IF
    1416             :          ELSE
    1417             :             ! Quadratic equation with max two roots.
    1418           0 :             DD = c*c - 4.0_dp*b*d
    1419           0 :             IF (DD .GT. 0.0_dp) THEN
    1420           0 :                nroots = 2
    1421           0 :                x(1) = (-c + SQRT(DD))/2.0_dp/b
    1422           0 :                x(2) = (-c - SQRT(DD))/2.0_dp/b
    1423           0 :             ELSE IF (DD .LT. 0.0_dp) THEN
    1424             :                nroots = 0
    1425             :             ELSE
    1426           0 :                nroots = 1
    1427           0 :                x(1) = -c/2.0_dp/b
    1428             :             END IF
    1429             :          END IF
    1430             :       ELSE
    1431             :          ! Cubic equation with max three roots
    1432             :          ! Calculate p and q
    1433           0 :          p = c/a - b*b/a/a/3.0_dp
    1434           0 :          q = (2.0_dp*b*b*b/a/a/a - 9.0_dp*b*c/a/a + 27.0_dp*d/a)/27.0_dp
    1435             : 
    1436             :          ! Calculate DD
    1437           0 :          DD = p*p*p/27.0_dp + q*q/4.0_dp
    1438             : 
    1439           0 :          IF (DD .LT. 0.0_dp) THEN
    1440             :             ! three real unequal roots -- use the trigonometric formulation
    1441           0 :             phi = ACOS(-q/2.0_dp/SQRT(ABS(p*p*p)/27.0_dp))
    1442           0 :             temp1 = 2.0_dp*SQRT(ABS(p)/3.0_dp)
    1443           0 :             y1 = temp1*COS(phi/3.0_dp)
    1444           0 :             y2 = -temp1*COS((phi + pi)/3.0_dp)
    1445           0 :             y3 = -temp1*COS((phi - pi)/3.0_dp)
    1446             :          ELSE
    1447             :             ! 1 real & 2 conjugate complex roots OR 3 real roots (some are equal)
    1448           0 :             temp1 = -q/2.0_dp + SQRT(DD)
    1449           0 :             temp2 = -q/2.0_dp - SQRT(DD)
    1450           0 :             u = ABS(temp1)**(1.0_dp/3.0_dp)
    1451           0 :             v = ABS(temp2)**(1.0_dp/3.0_dp)
    1452           0 :             IF (temp1 .LT. 0.0_dp) u = -u
    1453           0 :             IF (temp2 .LT. 0.0_dp) v = -v
    1454           0 :             y1 = u + v
    1455           0 :             y2r = -(u + v)/2.0_dp
    1456           0 :             y2i = (u - v)*SQRT(3.0_dp)/2.0_dp
    1457             :          END IF
    1458             : 
    1459             :          ! Final transformation
    1460           0 :          temp1 = b/a/3.0_dp
    1461           0 :          y1 = y1 - temp1
    1462           0 :          y2 = y2 - temp1
    1463           0 :          y3 = y3 - temp1
    1464           0 :          y2r = y2r - temp1
    1465             : 
    1466             :          ! Assign answers
    1467           0 :          IF (DD .LT. 0.0_dp) THEN
    1468           0 :             nroots = 3
    1469           0 :             x(1) = y1
    1470           0 :             x(2) = y2
    1471           0 :             x(3) = y3
    1472           0 :          ELSE IF (DD .EQ. 0.0_dp) THEN
    1473           0 :             nroots = 2
    1474           0 :             x(1) = y1
    1475           0 :             x(2) = y2r
    1476             :             !x(3) = cmplx(y2r,  0.)
    1477             :          ELSE
    1478           0 :             nroots = 1
    1479           0 :             x(1) = y1
    1480             :             !x(2) = cmplx(y2r, y2i)
    1481             :             !x(3) = cmplx(y2r,-y2i)
    1482             :          END IF
    1483             : 
    1484             :       END IF
    1485             : 
    1486             : !write(*,'(i2,a)') nroots, ' real root(s)'
    1487           0 :       nmins = 0
    1488           0 :       DO i = 1, nroots
    1489             :          ! maximum or minimum? use the derivative
    1490             :          ! 3*a*x**2+2*b*x+c
    1491           0 :          der = 3.0_dp*a*x(i)*x(i) + 2.0_dp*b*x(i) + c
    1492           0 :          IF (der .GT. 0.0_dp) THEN
    1493           0 :             nmins = nmins + 1
    1494           0 :             minima(nmins) = x(i)
    1495             : !write(*,'(a,i2,a,f10.5)') 'Minimum ', i, ', value: ', x(i)
    1496             :          END IF
    1497             :       END DO
    1498             : 
    1499             : !    CALL timestop(handle)
    1500             : 
    1501           0 :    END SUBROUTINE analytic_line_search
    1502             : 
    1503             : ! **************************************************************************************************
    1504             : !> \brief Diagonalizes diagonal blocks of a symmetric dbcsr matrix
    1505             : !>        and returs its eigenvectors
    1506             : !> \param matrix ...
    1507             : !> \param c ...
    1508             : !> \param e ...
    1509             : !> \par History
    1510             : !>       2011.07 created [Rustam Z Khaliullin]
    1511             : !> \author Rustam Z Khaliullin
    1512             : ! **************************************************************************************************
    1513           0 :    SUBROUTINE diagonalize_diagonal_blocks(matrix, c, e)
    1514             : 
    1515             :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix
    1516             :       TYPE(dbcsr_type), INTENT(OUT)                      :: c
    1517             :       TYPE(dbcsr_type), INTENT(OUT), OPTIONAL            :: e
    1518             : 
    1519             :       CHARACTER(len=*), PARAMETER :: routineN = 'diagonalize_diagonal_blocks'
    1520             : 
    1521             :       INTEGER                                            :: handle, iblock_col, iblock_row, &
    1522             :                                                             iblock_size, info, lwork, orbital
    1523             :       LOGICAL                                            :: block_needed, do_eigenvalues
    1524           0 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues, work
    1525           0 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: data_copy
    1526           0 :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: data_p, p_new_block
    1527             :       TYPE(dbcsr_iterator_type)                          :: iter
    1528             : 
    1529           0 :       CALL timeset(routineN, handle)
    1530             : 
    1531           0 :       IF (PRESENT(e)) THEN
    1532             :          do_eigenvalues = .TRUE.
    1533             :       ELSE
    1534           0 :          do_eigenvalues = .FALSE.
    1535             :       END IF
    1536             : 
    1537             :       ! create a matrix for eigenvectors
    1538           0 :       CALL dbcsr_work_create(c, work_mutable=.TRUE.)
    1539           0 :       IF (do_eigenvalues) &
    1540           0 :          CALL dbcsr_work_create(e, work_mutable=.TRUE.)
    1541             : 
    1542           0 :       CALL dbcsr_iterator_start(iter, matrix)
    1543             : 
    1544           0 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
    1545             : 
    1546           0 :          CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
    1547             : 
    1548           0 :          block_needed = .FALSE.
    1549           0 :          IF (iblock_row == iblock_col) block_needed = .TRUE.
    1550             : 
    1551           0 :          IF (block_needed) THEN
    1552             : 
    1553             :             ! Prepare data
    1554           0 :             ALLOCATE (eigenvalues(iblock_size))
    1555           0 :             ALLOCATE (data_copy(iblock_size, iblock_size))
    1556           0 :             data_copy(:, :) = data_p(:, :)
    1557             : 
    1558             :             ! Query the optimal workspace for dsyev
    1559           0 :             LWORK = -1
    1560           0 :             ALLOCATE (WORK(MAX(1, LWORK)))
    1561           0 :             CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
    1562           0 :             LWORK = INT(WORK(1))
    1563           0 :             DEALLOCATE (WORK)
    1564             : 
    1565             :             ! Allocate the workspace and solve the eigenproblem
    1566           0 :             ALLOCATE (WORK(MAX(1, LWORK)))
    1567           0 :             CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
    1568           0 :             IF (INFO .NE. 0) CPABORT("DSYEV failed")
    1569             : 
    1570             :             ! copy eigenvectors into a cp_dbcsr matrix
    1571           0 :             NULLIFY (p_new_block)
    1572           0 :             CALL dbcsr_reserve_block2d(c, iblock_row, iblock_col, p_new_block)
    1573           0 :             CPASSERT(ASSOCIATED(p_new_block))
    1574           0 :             p_new_block(:, :) = data_copy(:, :)
    1575             : 
    1576             :             ! if requested copy eigenvalues into a cp_dbcsr matrix
    1577           0 :             IF (do_eigenvalues) THEN
    1578           0 :                NULLIFY (p_new_block)
    1579           0 :                CALL dbcsr_reserve_block2d(e, iblock_row, iblock_col, p_new_block)
    1580           0 :                CPASSERT(ASSOCIATED(p_new_block))
    1581           0 :                p_new_block(:, :) = 0.0_dp
    1582           0 :                DO orbital = 1, iblock_size
    1583           0 :                   p_new_block(orbital, orbital) = eigenvalues(orbital)
    1584             :                END DO
    1585             :             END IF
    1586             : 
    1587           0 :             DEALLOCATE (WORK)
    1588           0 :             DEALLOCATE (data_copy)
    1589           0 :             DEALLOCATE (eigenvalues)
    1590             : 
    1591             :          END IF
    1592             : 
    1593             :       END DO
    1594             : 
    1595           0 :       CALL dbcsr_iterator_stop(iter)
    1596             : 
    1597           0 :       CALL dbcsr_finalize(c)
    1598           0 :       IF (do_eigenvalues) CALL dbcsr_finalize(e)
    1599             : 
    1600           0 :       CALL timestop(handle)
    1601             : 
    1602           0 :    END SUBROUTINE diagonalize_diagonal_blocks
    1603             : 
    1604             : ! **************************************************************************************************
    1605             : !> \brief Transforms a matrix M_out = tr(U1) * M_in * U2
    1606             : !> \param matrix ...
    1607             : !> \param u1 ...
    1608             : !> \param u2 ...
    1609             : !> \param eps_filter ...
    1610             : !> \par History
    1611             : !>       2011.10 created [Rustam Z Khaliullin]
    1612             : !> \author Rustam Z Khaliullin
    1613             : ! **************************************************************************************************
    1614           0 :    SUBROUTINE matrix_forward_transform(matrix, u1, u2, eps_filter)
    1615             : 
    1616             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
    1617             :       TYPE(dbcsr_type), INTENT(IN)                       :: u1, u2
    1618             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
    1619             : 
    1620             :       CHARACTER(len=*), PARAMETER :: routineN = 'matrix_forward_transform'
    1621             : 
    1622             :       INTEGER                                            :: handle
    1623             :       TYPE(dbcsr_type)                                   :: tmp
    1624             : 
    1625           0 :       CALL timeset(routineN, handle)
    1626             : 
    1627             :       CALL dbcsr_create(tmp, template=matrix, &
    1628           0 :                         matrix_type=dbcsr_type_no_symmetry)
    1629             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, u2, 0.0_dp, tmp, &
    1630           0 :                           filter_eps=eps_filter)
    1631             :       CALL dbcsr_multiply("T", "N", 1.0_dp, u1, tmp, 0.0_dp, matrix, &
    1632           0 :                           filter_eps=eps_filter)
    1633           0 :       CALL dbcsr_release(tmp)
    1634             : 
    1635           0 :       CALL timestop(handle)
    1636             : 
    1637           0 :    END SUBROUTINE matrix_forward_transform
    1638             : 
    1639             : ! **************************************************************************************************
    1640             : !> \brief Transforms a matrix M_out = U1 * M_in * tr(U2)
    1641             : !> \param matrix ...
    1642             : !> \param u1 ...
    1643             : !> \param u2 ...
    1644             : !> \param eps_filter ...
    1645             : !> \par History
    1646             : !>       2011.10 created [Rustam Z Khaliullin]
    1647             : !> \author Rustam Z Khaliullin
    1648             : ! **************************************************************************************************
    1649           0 :    SUBROUTINE matrix_backward_transform(matrix, u1, u2, eps_filter)
    1650             : 
    1651             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
    1652             :       TYPE(dbcsr_type), INTENT(IN)                       :: u1, u2
    1653             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
    1654             : 
    1655             :       CHARACTER(len=*), PARAMETER :: routineN = 'matrix_backward_transform'
    1656             : 
    1657             :       INTEGER                                            :: handle
    1658             :       TYPE(dbcsr_type)                                   :: tmp
    1659             : 
    1660           0 :       CALL timeset(routineN, handle)
    1661             : 
    1662             :       CALL dbcsr_create(tmp, template=matrix, &
    1663           0 :                         matrix_type=dbcsr_type_no_symmetry)
    1664             :       CALL dbcsr_multiply("N", "T", 1.0_dp, matrix, u2, 0.0_dp, tmp, &
    1665           0 :                           filter_eps=eps_filter)
    1666             :       CALL dbcsr_multiply("N", "N", 1.0_dp, u1, tmp, 0.0_dp, matrix, &
    1667           0 :                           filter_eps=eps_filter)
    1668           0 :       CALL dbcsr_release(tmp)
    1669             : 
    1670           0 :       CALL timestop(handle)
    1671             : 
    1672           0 :    END SUBROUTINE matrix_backward_transform
    1673             : 
    1674             : !! **************************************************************************************************
    1675             : !!> \brief Transforms to a representation in which diagonal blocks
    1676             : !!>        of qq and pp matrices are diagonal. This can improve convergence
    1677             : !!>        of PCG
    1678             : !!> \par History
    1679             : !!>       2011.07 created [Rustam Z Khaliullin]
    1680             : !!> \author Rustam Z Khaliullin
    1681             : !! **************************************************************************************************
    1682             : !  SUBROUTINE transform_matrices_to_blk_diag(matrix_pp,matrix_qq,matrix_qp,&
    1683             : !    matrix_pq,eps_filter)
    1684             : !
    1685             : !    TYPE(dbcsr_type), INTENT(INOUT)       :: matrix_pp, matrix_qq,&
    1686             : !                                                matrix_qp, matrix_pq
    1687             : !    REAL(KIND=dp), INTENT(IN)                :: eps_filter
    1688             : !
    1689             : !    CHARACTER(len=*), PARAMETER :: routineN = 'transform_matrices_to_blk_diag',&
    1690             : !      routineP = moduleN//':'//routineN
    1691             : !
    1692             : !    TYPE(dbcsr_type)                      :: tmp_pp, tmp_qq,&
    1693             : !                                                tmp_qp, tmp_pq,&
    1694             : !                                                blk, blk2
    1695             : !    INTEGER                                  :: handle
    1696             : !
    1697             : !    CALL timeset(routineN,handle)
    1698             : !
    1699             : !    ! find a better basis by diagonalizing diagonal blocks
    1700             : !    ! first pp
    1701             : !    CALL dbcsr_init(blk)
    1702             : !    CALL dbcsr_create(blk,template=matrix_pp)
    1703             : !    CALL diagonalize_diagonal_blocks(matrix_pp,blk)
    1704             : !
    1705             : !    ! convert matrices to the new basis
    1706             : !    CALL dbcsr_init(tmp_pp)
    1707             : !    CALL dbcsr_create(tmp_pp,template=matrix_pp)
    1708             : !    CALL dbcsr_multiply("N","N",1.0_dp,matrix_pp,blk,0.0_dp,tmp_pp,&
    1709             : !               filter_eps=eps_filter)
    1710             : !    CALL dbcsr_multiply("T","N",1.0_dp,blk,tmp_pp,0.0_dp,matrix_pp,&
    1711             : !               filter_eps=eps_filter)
    1712             : !    CALL dbcsr_release(tmp_pp)
    1713             : !
    1714             : !    ! now qq
    1715             : !    CALL dbcsr_init(blk2)
    1716             : !    CALL dbcsr_create(blk2,template=matrix_qq)
    1717             : !    CALL diagonalize_diagonal_blocks(matrix_qq,blk2)
    1718             : !
    1719             : !    CALL dbcsr_init(tmp_qq)
    1720             : !    CALL dbcsr_create(tmp_qq,template=matrix_qq)
    1721             : !    CALL dbcsr_multiply("N","N",1.0_dp,matrix_qq,blk2,0.0_dp,tmp_qq,&
    1722             : !               filter_eps=eps_filter)
    1723             : !    CALL dbcsr_multiply("T","N",1.0_dp,blk2,tmp_qq,0.0_dp,matrix_qq,&
    1724             : !               filter_eps=eps_filter)
    1725             : !    CALL dbcsr_release(tmp_qq)
    1726             : !
    1727             : !    ! transform pq
    1728             : !    CALL dbcsr_init(tmp_pq)
    1729             : !    CALL dbcsr_create(tmp_pq,template=matrix_pq)
    1730             : !    CALL dbcsr_multiply("T","N",1.0_dp,blk,matrix_pq,0.0_dp,tmp_pq,&
    1731             : !               filter_eps=eps_filter)
    1732             : !    CALL dbcsr_multiply("N","N",1.0_dp,tmp_pq,blk2,0.0_dp,matrix_pq,&
    1733             : !               filter_eps=eps_filter)
    1734             : !    CALL dbcsr_release(tmp_pq)
    1735             : !
    1736             : !    ! transform qp
    1737             : !    CALL dbcsr_init(tmp_qp)
    1738             : !    CALL dbcsr_create(tmp_qp,template=matrix_qp)
    1739             : !    CALL dbcsr_multiply("N","N",1.0_dp,matrix_qp,blk,0.0_dp,tmp_qp,&
    1740             : !               filter_eps=eps_filter)
    1741             : !    CALL dbcsr_multiply("T","N",1.0_dp,blk2,tmp_qp,0.0_dp,matrix_qp,&
    1742             : !               filter_eps=eps_filter)
    1743             : !    CALL dbcsr_release(tmp_qp)
    1744             : !
    1745             : !    CALL dbcsr_release(blk2)
    1746             : !    CALL dbcsr_release(blk)
    1747             : !
    1748             : !    CALL timestop(handle)
    1749             : !
    1750             : !  END SUBROUTINE transform_matrices_to_blk_diag
    1751             : 
    1752             : ! **************************************************************************************************
    1753             : !> \brief computes oo, ov, vo, and vv blocks of the ks matrix
    1754             : !> \par History
    1755             : !>       2011.06 created [Rustam Z Khaliullin]
    1756             : !> \author Rustam Z Khaliullin
    1757             : ! **************************************************************************************************
    1758             : !  SUBROUTINE ct_step_env_execute(env)
    1759             : !
    1760             : !    TYPE(ct_step_env_type)                      :: env
    1761             : !
    1762             : !    CHARACTER(len=*), PARAMETER :: routineN = 'ct_step_env_execute', &
    1763             : !      routineP = moduleN//':'//routineN
    1764             : !
    1765             : !    INTEGER                                  :: handle
    1766             : !
    1767             : !    CALL timeset(routineN,handle)
    1768             : !
    1769             : !
    1770             : !    CALL timestop(handle)
    1771             : !
    1772             : !  END SUBROUTINE ct_step_env_execute
    1773             : 
    1774             : END MODULE ct_methods
    1775             : 

Generated by: LCOV version 1.15