LCOV - code coverage report
Current view: top level - src - ct_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 0 554 0.0 %
Date: 2024-03-28 07:31:50 Functions: 0 8 0.0 %

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

Generated by: LCOV version 1.15