LCOV - code coverage report
Current view: top level - src - mao_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 279 365 76.4 %
Date: 2025-02-18 08:24:35 Functions: 9 11 81.8 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Calculate MAO's and analyze wavefunctions
      10             : !> \par History
      11             : !>      03.2016 created [JGH]
      12             : !>      12.2016 split into four modules [JGH]
      13             : !> \author JGH
      14             : ! **************************************************************************************************
      15             : MODULE mao_methods
      16             :    USE atomic_kind_types,               ONLY: get_atomic_kind
      17             :    USE basis_set_container_types,       ONLY: add_basis_set_to_container
      18             :    USE basis_set_types,                 ONLY: create_primitive_basis_set,&
      19             :                                               get_gto_basis_set,&
      20             :                                               gto_basis_set_p_type,&
      21             :                                               gto_basis_set_type,&
      22             :                                               write_gto_basis_set
      23             :    USE cp_control_types,                ONLY: dft_control_type
      24             :    USE cp_dbcsr_api,                    ONLY: &
      25             :         dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_dot, dbcsr_get_block_p, &
      26             :         dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      27             :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
      28             :         dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_set, dbcsr_type, &
      29             :         dbcsr_type_no_symmetry
      30             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      31             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      32             :                                               cp_dbcsr_plus_fm_fm_t,&
      33             :                                               dbcsr_allocate_matrix_set
      34             :    USE cp_fm_diag,                      ONLY: cp_fm_geeig
      35             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      36             :                                               cp_fm_struct_release,&
      37             :                                               cp_fm_struct_type
      38             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      39             :                                               cp_fm_release,&
      40             :                                               cp_fm_type
      41             :    USE input_constants,                 ONLY: mao_basis_ext,&
      42             :                                               mao_basis_orb,&
      43             :                                               mao_basis_prim
      44             :    USE iterate_matrix,                  ONLY: invert_Hotelling
      45             :    USE kinds,                           ONLY: dp
      46             :    USE kpoint_methods,                  ONLY: rskp_transform
      47             :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      48             :                                               kpoint_type
      49             :    USE message_passing,                 ONLY: mp_comm_type,&
      50             :                                               mp_para_env_type
      51             :    USE particle_types,                  ONLY: particle_type
      52             :    USE qs_environment_types,            ONLY: get_qs_env,&
      53             :                                               qs_environment_type
      54             :    USE qs_interactions,                 ONLY: init_interaction_radii_orb_basis
      55             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      56             :                                               qs_kind_type
      57             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      58             : #include "./base/base_uses.f90"
      59             : 
      60             :    IMPLICIT NONE
      61             :    PRIVATE
      62             : 
      63             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_methods'
      64             : 
      65             :    TYPE mblocks
      66             :       INTEGER                                        :: n = -1, ma = -1
      67             :       REAL(KIND=dp), DIMENSION(:, :), POINTER        :: mat => NULL()
      68             :       REAL(KIND=dp), DIMENSION(:), POINTER           :: eig => NULL()
      69             :    END TYPE mblocks
      70             : 
      71             :    PUBLIC :: mao_initialization, mao_function, mao_function_gradient, mao_orthogonalization, &
      72             :              mao_project_gradient, mao_scalar_product, mao_build_q, mao_basis_analysis, &
      73             :              mao_reference_basis, calculate_p_gamma
      74             : 
      75             : ! **************************************************************************************************
      76             : 
      77             : CONTAINS
      78             : 
      79             : ! **************************************************************************************************
      80             : !> \brief ...
      81             : !> \param mao_coef ...
      82             : !> \param pmat ...
      83             : !> \param smat ...
      84             : !> \param eps1 ...
      85             : !> \param iolevel ...
      86             : !> \param iw ...
      87             : ! **************************************************************************************************
      88          16 :    SUBROUTINE mao_initialization(mao_coef, pmat, smat, eps1, iolevel, iw)
      89             :       TYPE(dbcsr_type)                                   :: mao_coef, pmat, smat
      90             :       REAL(KIND=dp), INTENT(IN)                          :: eps1
      91             :       INTEGER, INTENT(IN)                                :: iolevel, iw
      92             : 
      93             :       INTEGER                                            :: i, iatom, info, jatom, lwork, m, n, nblk
      94          16 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_sizes, mao_blk, row_blk, &
      95          16 :                                                             row_blk_sizes
      96             :       LOGICAL                                            :: found
      97          16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: w, work
      98          16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: amat, bmat
      99          16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cblock, pblock, sblock
     100             :       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
     101             :       TYPE(dbcsr_iterator_type)                          :: dbcsr_iter
     102          16 :       TYPE(mblocks), ALLOCATABLE, DIMENSION(:)           :: mbl
     103             :       TYPE(mp_comm_type)                                 :: group
     104             : 
     105          16 :       CALL dbcsr_get_info(mao_coef, nblkrows_total=nblk)
     106         126 :       ALLOCATE (mbl(nblk))
     107          94 :       DO i = 1, nblk
     108          94 :          NULLIFY (mbl(i)%mat, mbl(i)%eig)
     109             :       END DO
     110             : 
     111          16 :       CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
     112          55 :       DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
     113          39 :          CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
     114          39 :          CPASSERT(iatom == jatom)
     115          39 :          m = SIZE(cblock, 2)
     116          39 :          NULLIFY (pblock, sblock)
     117          39 :          CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pblock, found=found)
     118          39 :          CPASSERT(found)
     119          39 :          CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
     120          39 :          CPASSERT(found)
     121          39 :          n = SIZE(sblock, 1)
     122          39 :          lwork = MAX(n*n, 100)
     123         390 :          ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork))
     124       20925 :          amat(1:n, 1:n) = pblock(1:n, 1:n)
     125       20925 :          bmat(1:n, 1:n) = sblock(1:n, 1:n)
     126          39 :          info = 0
     127          39 :          CALL dsygv(1, "V", "U", n, amat, n, bmat, n, w, work, lwork, info)
     128          39 :          CPASSERT(info == 0)
     129         234 :          ALLOCATE (mbl(iatom)%mat(n, n), mbl(iatom)%eig(n))
     130          39 :          mbl(iatom)%n = n
     131          39 :          mbl(iatom)%ma = m
     132         797 :          DO i = 1, n
     133         758 :             mbl(iatom)%eig(i) = w(n - i + 1)
     134       20925 :             mbl(iatom)%mat(1:n, i) = amat(1:n, n - i + 1)
     135             :          END DO
     136        2105 :          cblock(1:n, 1:m) = amat(1:n, n:n - m + 1:-1)
     137         172 :          DEALLOCATE (amat, bmat, w, work)
     138             :       END DO
     139          16 :       CALL dbcsr_iterator_stop(dbcsr_iter)
     140             : 
     141          16 :       IF (eps1 < 10.0_dp) THEN
     142           0 :          CALL dbcsr_get_info(mao_coef, row_blk_size=row_blk_sizes, group=group)
     143           0 :          ALLOCATE (row_blk(nblk), mao_blk(nblk))
     144           0 :          mao_blk = 0
     145           0 :          row_blk = row_blk_sizes
     146           0 :          DO iatom = 1, nblk
     147           0 :             IF (ASSOCIATED(mbl(iatom)%mat)) THEN
     148           0 :                n = mbl(iatom)%n
     149           0 :                m = 0
     150           0 :                DO i = 1, n
     151           0 :                   IF (mbl(iatom)%eig(i) < eps1) EXIT
     152           0 :                   m = i
     153             :                END DO
     154           0 :                m = MAX(m, mbl(iatom)%ma)
     155           0 :                mbl(iatom)%ma = m
     156           0 :                mao_blk(iatom) = m
     157             :             END IF
     158             :          END DO
     159           0 :          CALL group%sum(mao_blk)
     160           0 :          CALL dbcsr_get_info(mao_coef, distribution=dbcsr_dist)
     161           0 :          CALL dbcsr_release(mao_coef)
     162             :          CALL dbcsr_create(mao_coef, name="MAO_COEF", dist=dbcsr_dist, &
     163             :                            matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_blk, &
     164           0 :                            col_blk_size=mao_blk, nze=0)
     165           0 :          CALL dbcsr_reserve_diag_blocks(matrix=mao_coef)
     166           0 :          DEALLOCATE (mao_blk, row_blk)
     167             :          !
     168           0 :          CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
     169           0 :          DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
     170           0 :             CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
     171           0 :             CPASSERT(iatom == jatom)
     172           0 :             n = SIZE(cblock, 1)
     173           0 :             m = SIZE(cblock, 2)
     174           0 :             CPASSERT(n == mbl(iatom)%n .AND. m == mbl(iatom)%ma)
     175           0 :             cblock(1:n, 1:m) = mbl(iatom)%mat(1:n, 1:m)
     176             :          END DO
     177           0 :          CALL dbcsr_iterator_stop(dbcsr_iter)
     178             :          !
     179             :       END IF
     180             : 
     181          16 :       IF (iolevel > 2) THEN
     182             :          CALL dbcsr_get_info(mao_coef, col_blk_size=col_blk_sizes, &
     183          12 :                              row_blk_size=row_blk_sizes, group=group)
     184          66 :          DO iatom = 1, nblk
     185          54 :             n = row_blk_sizes(iatom)
     186          54 :             m = col_blk_sizes(iatom)
     187         162 :             ALLOCATE (w(n))
     188         978 :             w(1:n) = 0._dp
     189          54 :             IF (ASSOCIATED(mbl(iatom)%mat)) THEN
     190         489 :                w(1:n) = mbl(iatom)%eig(1:n)
     191             :             END IF
     192          54 :             CALL group%sum(w)
     193          54 :             IF (iw > 0) THEN
     194          27 :                WRITE (iw, '(A,i2,20F8.4)', ADVANCE="NO") " Spectrum/Gap  ", iatom, w(1:m)
     195          27 :                WRITE (iw, '(A,F8.4)') " || ", w(m + 1)
     196             :             END IF
     197          66 :             DEALLOCATE (w)
     198             :          END DO
     199             :       END IF
     200             : 
     201          16 :       CALL mao_orthogonalization(mao_coef, smat)
     202             : 
     203          94 :       DO i = 1, nblk
     204          78 :          IF (ASSOCIATED(mbl(i)%mat)) THEN
     205          39 :             DEALLOCATE (mbl(i)%mat)
     206             :          END IF
     207          94 :          IF (ASSOCIATED(mbl(i)%eig)) THEN
     208          39 :             DEALLOCATE (mbl(i)%eig)
     209             :          END IF
     210             :       END DO
     211          16 :       DEALLOCATE (mbl)
     212             : 
     213          48 :    END SUBROUTINE mao_initialization
     214             : 
     215             : ! **************************************************************************************************
     216             : !> \brief ...
     217             : !> \param mao_coef ...
     218             : !> \param fval ...
     219             : !> \param qmat ...
     220             : !> \param smat ...
     221             : !> \param binv ...
     222             : !> \param reuse ...
     223             : ! **************************************************************************************************
     224         318 :    SUBROUTINE mao_function(mao_coef, fval, qmat, smat, binv, reuse)
     225             :       TYPE(dbcsr_type)                                   :: mao_coef
     226             :       REAL(KIND=dp), INTENT(OUT)                         :: fval
     227             :       TYPE(dbcsr_type)                                   :: qmat, smat, binv
     228             :       LOGICAL, INTENT(IN)                                :: reuse
     229             : 
     230             :       REAL(KIND=dp)                                      :: convergence, threshold
     231             :       TYPE(dbcsr_type)                                   :: bmat, scmat, tmat
     232             : 
     233         318 :       threshold = 1.e-8_dp
     234         318 :       convergence = 1.e-6_dp
     235             :       ! temp matrices
     236         318 :       CALL dbcsr_create(scmat, template=mao_coef)
     237         318 :       CALL dbcsr_create(bmat, template=binv)
     238         318 :       CALL dbcsr_create(tmat, template=qmat)
     239             :       ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
     240         318 :       CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
     241         318 :       CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
     242             :       ! calculate inverse of B
     243             :       CALL invert_Hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
     244         318 :                             norm_convergence=convergence, silent=.TRUE.)
     245             :       ! calculate Binv*C and T=C(T)*Binv*C
     246         318 :       CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
     247         318 :       CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
     248             :       ! function = Tr(Q*T)
     249         318 :       CALL dbcsr_dot(qmat, tmat, fval)
     250             :       ! free temp matrices
     251         318 :       CALL dbcsr_release(scmat)
     252         318 :       CALL dbcsr_release(bmat)
     253         318 :       CALL dbcsr_release(tmat)
     254             : 
     255         318 :    END SUBROUTINE mao_function
     256             : 
     257             : ! **************************************************************************************************
     258             : !> \brief ...
     259             : !> \param mao_coef ...
     260             : !> \param fval ...
     261             : !> \param mao_grad ...
     262             : !> \param qmat ...
     263             : !> \param smat ...
     264             : !> \param binv ...
     265             : !> \param reuse ...
     266             : ! **************************************************************************************************
     267         166 :    SUBROUTINE mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse)
     268             :       TYPE(dbcsr_type)                                   :: mao_coef
     269             :       REAL(KIND=dp), INTENT(OUT)                         :: fval
     270             :       TYPE(dbcsr_type)                                   :: mao_grad, qmat, smat, binv
     271             :       LOGICAL, INTENT(IN)                                :: reuse
     272             : 
     273             :       REAL(KIND=dp)                                      :: convergence, threshold
     274             :       TYPE(dbcsr_type)                                   :: bmat, scmat, t2mat, tmat
     275             : 
     276         166 :       threshold = 1.e-8_dp
     277         166 :       convergence = 1.e-6_dp
     278             :       ! temp matrices
     279         166 :       CALL dbcsr_create(scmat, template=mao_coef)
     280         166 :       CALL dbcsr_create(bmat, template=binv)
     281         166 :       CALL dbcsr_create(tmat, template=qmat)
     282         166 :       CALL dbcsr_create(t2mat, template=scmat)
     283             :       ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap
     284         166 :       CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat)
     285         166 :       CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat)
     286             :       ! calculate inverse of B
     287             :       CALL invert_Hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, &
     288         166 :                             norm_convergence=convergence, silent=.TRUE.)
     289             :       ! calculate R=C*Binv and T=C*Binv*C(T)=R*C(T)
     290         166 :       CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat)
     291         166 :       CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat)
     292             :       ! function = Tr(Q*T)
     293         166 :       CALL dbcsr_dot(qmat, tmat, fval)
     294             :       ! Gradient part 1: g = 2*Q*C*Binv = 2*Q*R
     295             :       CALL dbcsr_multiply("N", "N", 2.0_dp, qmat, scmat, 0.0_dp, mao_grad, &
     296         166 :                           retain_sparsity=.TRUE.)
     297             :       ! Gradient part 2: g = -2*S*T*X; X = Q*R
     298         166 :       CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, scmat, 0.0_dp, t2mat)
     299         166 :       CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, t2mat, 0.0_dp, scmat)
     300             :       CALL dbcsr_multiply("N", "N", -2.0_dp, smat, scmat, 1.0_dp, mao_grad, &
     301         166 :                           retain_sparsity=.TRUE.)
     302             :       ! free temp matrices
     303         166 :       CALL dbcsr_release(scmat)
     304         166 :       CALL dbcsr_release(bmat)
     305         166 :       CALL dbcsr_release(tmat)
     306         166 :       CALL dbcsr_release(t2mat)
     307             : 
     308         166 :       CALL mao_project_gradient(mao_coef, mao_grad, smat)
     309             : 
     310         166 :    END SUBROUTINE mao_function_gradient
     311             : 
     312             : ! **************************************************************************************************
     313             : !> \brief ...
     314             : !> \param mao_coef ...
     315             : !> \param smat ...
     316             : ! **************************************************************************************************
     317         484 :    SUBROUTINE mao_orthogonalization(mao_coef, smat)
     318             :       TYPE(dbcsr_type)                                   :: mao_coef, smat
     319             : 
     320             :       INTEGER                                            :: i, iatom, info, jatom, lwork, m, n
     321             :       LOGICAL                                            :: found
     322         484 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: w, work
     323         484 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: amat, bmat
     324         484 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cblock, sblock
     325             :       TYPE(dbcsr_iterator_type)                          :: dbcsr_iter
     326             : 
     327         484 :       CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
     328        1747 :       DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
     329        1263 :          CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
     330        1263 :          CPASSERT(iatom == jatom)
     331        1263 :          m = SIZE(cblock, 2)
     332        1263 :          n = SIZE(cblock, 1)
     333        1263 :          NULLIFY (sblock)
     334        1263 :          CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
     335        1263 :          CPASSERT(found)
     336        1263 :          lwork = MAX(n*n, 100)
     337       13893 :          ALLOCATE (amat(n, m), bmat(m, m), w(m), work(lwork))
     338     2667347 :          amat(1:n, 1:m) = MATMUL(sblock(1:n, 1:n), cblock(1:n, 1:m))
     339      278171 :          bmat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), amat(1:n, 1:m))
     340        1263 :          info = 0
     341        1263 :          CALL dsyev("V", "U", m, bmat, m, w, work, lwork, info)
     342        1263 :          CPASSERT(info == 0)
     343        3789 :          CPASSERT(ALL(w > 0.0_dp))
     344        3789 :          w = 1.0_dp/SQRT(w)
     345        3789 :          DO i = 1, m
     346       11367 :             amat(1:m, i) = bmat(1:m, i)*w(i)
     347             :          END DO
     348       59361 :          bmat(1:m, 1:m) = MATMUL(amat(1:m, 1:m), TRANSPOSE(bmat(1:m, 1:m)))
     349      422583 :          cblock(1:n, 1:m) = MATMUL(cblock(1:n, 1:m), bmat(1:m, 1:m))
     350        3789 :          DEALLOCATE (amat, bmat, w, work)
     351             :       END DO
     352         484 :       CALL dbcsr_iterator_stop(dbcsr_iter)
     353             : 
     354         968 :    END SUBROUTINE mao_orthogonalization
     355             : 
     356             : ! **************************************************************************************************
     357             : !> \brief ...
     358             : !> \param mao_coef ...
     359             : !> \param mao_grad ...
     360             : !> \param smat ...
     361             : ! **************************************************************************************************
     362         166 :    SUBROUTINE mao_project_gradient(mao_coef, mao_grad, smat)
     363             :       TYPE(dbcsr_type)                                   :: mao_coef, mao_grad, smat
     364             : 
     365             :       INTEGER                                            :: iatom, jatom, m, n
     366             :       LOGICAL                                            :: found
     367         166 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: amat
     368         166 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cblock, gblock, sblock
     369             :       TYPE(dbcsr_iterator_type)                          :: dbcsr_iter
     370             : 
     371         166 :       CALL dbcsr_iterator_start(dbcsr_iter, mao_coef)
     372         598 :       DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
     373         432 :          CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock)
     374         432 :          CPASSERT(iatom == jatom)
     375         432 :          m = SIZE(cblock, 2)
     376         432 :          n = SIZE(cblock, 1)
     377         432 :          NULLIFY (sblock)
     378         432 :          CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found)
     379         432 :          CPASSERT(found)
     380         432 :          NULLIFY (gblock)
     381         432 :          CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=jatom, block=gblock, found=found)
     382         432 :          CPASSERT(found)
     383        1728 :          ALLOCATE (amat(m, m))
     384     1950960 :          amat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), MATMUL(sblock(1:n, 1:n), gblock(1:n, 1:m)))
     385      146304 :          gblock(1:n, 1:m) = gblock(1:n, 1:m) - MATMUL(cblock(1:n, 1:m), amat(1:m, 1:m))
     386        1728 :          DEALLOCATE (amat)
     387             :       END DO
     388         166 :       CALL dbcsr_iterator_stop(dbcsr_iter)
     389             : 
     390         332 :    END SUBROUTINE mao_project_gradient
     391             : 
     392             : ! **************************************************************************************************
     393             : !> \brief ...
     394             : !> \param fmat1 ...
     395             : !> \param fmat2 ...
     396             : !> \return ...
     397             : ! **************************************************************************************************
     398         332 :    FUNCTION mao_scalar_product(fmat1, fmat2) RESULT(spro)
     399             :       TYPE(dbcsr_type)                                   :: fmat1, fmat2
     400             :       REAL(KIND=dp)                                      :: spro
     401             : 
     402             :       INTEGER                                            :: iatom, jatom, m, n
     403             :       LOGICAL                                            :: found
     404         166 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ablock, bblock
     405             :       TYPE(dbcsr_iterator_type)                          :: dbcsr_iter
     406             :       TYPE(mp_comm_type)                                 :: group
     407             : 
     408         166 :       spro = 0.0_dp
     409             : 
     410         166 :       CALL dbcsr_iterator_start(dbcsr_iter, fmat1)
     411         598 :       DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
     412         432 :          CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, ablock)
     413         432 :          CPASSERT(iatom == jatom)
     414         432 :          m = SIZE(ablock, 2)
     415         432 :          n = SIZE(ablock, 1)
     416         432 :          CALL dbcsr_get_block_p(matrix=fmat2, row=iatom, col=jatom, block=bblock, found=found)
     417         432 :          CPASSERT(found)
     418       26950 :          spro = spro + SUM(ablock(1:n, 1:m)*bblock(1:n, 1:m))
     419             :       END DO
     420         166 :       CALL dbcsr_iterator_stop(dbcsr_iter)
     421             : 
     422         166 :       CALL dbcsr_get_info(fmat1, group=group)
     423         166 :       CALL group%sum(spro)
     424             : 
     425         166 :    END FUNCTION mao_scalar_product
     426             : 
     427             : ! **************************************************************************************************
     428             : !> \brief Calculate the density matrix at the Gamma point
     429             : !> \param pmat ...
     430             : !> \param ksmat ...
     431             : !> \param smat ...
     432             : !> \param kpoints      Kpoint environment
     433             : !> \param nmos         Number of occupied orbitals
     434             : !> \param occ          Maximum occupation per orbital
     435             : !> \par History
     436             : !>      04.2016 created [JGH]
     437             : ! **************************************************************************************************
     438           0 :    SUBROUTINE calculate_p_gamma(pmat, ksmat, smat, kpoints, nmos, occ)
     439             : 
     440             :       TYPE(dbcsr_type)                                   :: pmat, ksmat, smat
     441             :       TYPE(kpoint_type), POINTER                         :: kpoints
     442             :       INTEGER, INTENT(IN)                                :: nmos
     443             :       REAL(KIND=dp), INTENT(IN)                          :: occ
     444             : 
     445             :       INTEGER                                            :: norb
     446             :       REAL(KIND=dp)                                      :: de
     447             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues
     448             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     449             :       TYPE(cp_fm_type)                                   :: fmksmat, fmsmat, fmvec, fmwork
     450             :       TYPE(dbcsr_type)                                   :: tempmat
     451             : 
     452             :       ! FM matrices
     453             : 
     454           0 :       CALL dbcsr_get_info(smat, nfullrows_total=norb)
     455             :       CALL cp_fm_struct_create(fmstruct=matrix_struct, context=kpoints%blacs_env_all, &
     456           0 :                                nrow_global=norb, ncol_global=norb)
     457           0 :       CALL cp_fm_create(fmksmat, matrix_struct)
     458           0 :       CALL cp_fm_create(fmsmat, matrix_struct)
     459           0 :       CALL cp_fm_create(fmvec, matrix_struct)
     460           0 :       CALL cp_fm_create(fmwork, matrix_struct)
     461           0 :       ALLOCATE (eigenvalues(norb))
     462             : 
     463             :       ! DBCSR matrix
     464           0 :       CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry)
     465             : 
     466             :       ! transfer to FM
     467           0 :       CALL dbcsr_desymmetrize(smat, tempmat)
     468           0 :       CALL copy_dbcsr_to_fm(tempmat, fmsmat)
     469           0 :       CALL dbcsr_desymmetrize(ksmat, tempmat)
     470           0 :       CALL copy_dbcsr_to_fm(tempmat, fmksmat)
     471             : 
     472             :       ! diagonalize
     473           0 :       CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork)
     474           0 :       de = eigenvalues(nmos + 1) - eigenvalues(nmos)
     475           0 :       IF (de < 0.001_dp) THEN
     476             :          CALL cp_warn(__LOCATION__, "MAO: No band gap at "// &
     477           0 :                       "Gamma point. MAO analysis not reliable.")
     478             :       END IF
     479             :       ! density matrix
     480           0 :       CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pmat, matrix_v=fmvec, ncol=nmos, alpha=occ)
     481             : 
     482           0 :       DEALLOCATE (eigenvalues)
     483           0 :       CALL dbcsr_release(tempmat)
     484           0 :       CALL cp_fm_release(fmksmat)
     485           0 :       CALL cp_fm_release(fmsmat)
     486           0 :       CALL cp_fm_release(fmvec)
     487           0 :       CALL cp_fm_release(fmwork)
     488           0 :       CALL cp_fm_struct_release(matrix_struct)
     489             : 
     490           0 :    END SUBROUTINE calculate_p_gamma
     491             : 
     492             : ! **************************************************************************************************
     493             : !> \brief Define the MAO reference basis set
     494             : !> \param qs_env ...
     495             : !> \param mao_basis ...
     496             : !> \param mao_basis_set_list ...
     497             : !> \param orb_basis_set_list ...
     498             : !> \param iunit ...
     499             : !> \param print_basis ...
     500             : !> \par History
     501             : !>      07.2016 created [JGH]
     502             : ! **************************************************************************************************
     503          10 :    SUBROUTINE mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, &
     504             :                                   iunit, print_basis)
     505             : 
     506             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     507             :       INTEGER, INTENT(IN)                                :: mao_basis
     508             :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: mao_basis_set_list, orb_basis_set_list
     509             :       INTEGER, INTENT(IN), OPTIONAL                      :: iunit
     510             :       LOGICAL, INTENT(IN), OPTIONAL                      :: print_basis
     511             : 
     512             :       INTEGER                                            :: ikind, nbas, nkind, unit_nr
     513             :       REAL(KIND=dp)                                      :: eps_pgf_orb
     514             :       TYPE(dft_control_type), POINTER                    :: dft_control
     515             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set, pbasis
     516          10 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     517             :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     518             : 
     519             :       ! Reference basis set
     520           0 :       CPASSERT(.NOT. ASSOCIATED(mao_basis_set_list))
     521          10 :       CPASSERT(.NOT. ASSOCIATED(orb_basis_set_list))
     522             : 
     523             :       ! options
     524          10 :       IF (PRESENT(iunit)) THEN
     525          10 :          unit_nr = iunit
     526             :       ELSE
     527           0 :          unit_nr = -1
     528             :       END IF
     529             : 
     530          10 :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
     531          10 :       nkind = SIZE(qs_kind_set)
     532          90 :       ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind))
     533          30 :       DO ikind = 1, nkind
     534          20 :          NULLIFY (mao_basis_set_list(ikind)%gto_basis_set)
     535          30 :          NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
     536             :       END DO
     537             :       !
     538          30 :       DO ikind = 1, nkind
     539          20 :          qs_kind => qs_kind_set(ikind)
     540          20 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
     541          30 :          IF (ASSOCIATED(basis_set)) orb_basis_set_list(ikind)%gto_basis_set => basis_set
     542             :       END DO
     543             :       !
     544          12 :       SELECT CASE (mao_basis)
     545             :       CASE (mao_basis_orb)
     546           6 :          DO ikind = 1, nkind
     547           4 :             qs_kind => qs_kind_set(ikind)
     548           4 :             CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
     549           6 :             IF (ASSOCIATED(basis_set)) mao_basis_set_list(ikind)%gto_basis_set => basis_set
     550             :          END DO
     551             :       CASE (mao_basis_prim)
     552           6 :          DO ikind = 1, nkind
     553           4 :             qs_kind => qs_kind_set(ikind)
     554           4 :             CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
     555           4 :             NULLIFY (pbasis)
     556           6 :             IF (ASSOCIATED(basis_set)) THEN
     557           4 :                CALL create_primitive_basis_set(basis_set, pbasis)
     558           4 :                CALL get_qs_env(qs_env, dft_control=dft_control)
     559           4 :                eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
     560           4 :                CALL init_interaction_radii_orb_basis(pbasis, eps_pgf_orb)
     561           4 :                pbasis%kind_radius = basis_set%kind_radius
     562           4 :                mao_basis_set_list(ikind)%gto_basis_set => pbasis
     563           4 :                CALL add_basis_set_to_container(qs_kind%basis_sets, pbasis, "MAO")
     564             :             END IF
     565             :          END DO
     566             :       CASE (mao_basis_ext)
     567          18 :          DO ikind = 1, nkind
     568          12 :             qs_kind => qs_kind_set(ikind)
     569          12 :             CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="MAO")
     570          18 :             IF (ASSOCIATED(basis_set)) THEN
     571          12 :                basis_set%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius
     572          12 :                mao_basis_set_list(ikind)%gto_basis_set => basis_set
     573             :             END IF
     574             :          END DO
     575             :       CASE DEFAULT
     576          10 :          CPABORT("Unknown option for MAO basis")
     577             :       END SELECT
     578          10 :       IF (unit_nr > 0) THEN
     579          15 :          DO ikind = 1, nkind
     580          15 :             IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN
     581             :                WRITE (UNIT=unit_nr, FMT="(T2,A,I4)") &
     582           0 :                   "WARNING: No MAO basis set associated with Kind ", ikind
     583             :             ELSE
     584          10 :                nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf
     585             :                WRITE (UNIT=unit_nr, FMT="(T2,A,I4,T56,A,I10)") &
     586          10 :                   "MAO basis set Kind ", ikind, " Number of BSF:", nbas
     587             :             END IF
     588             :          END DO
     589             :       END IF
     590             : 
     591          10 :       IF (PRESENT(print_basis)) THEN
     592          10 :          IF (print_basis) THEN
     593           0 :             DO ikind = 1, nkind
     594           0 :                basis_set => mao_basis_set_list(ikind)%gto_basis_set
     595           0 :                IF (ASSOCIATED(basis_set)) CALL write_gto_basis_set(basis_set, unit_nr, "MAO REFERENCE BASIS")
     596             :             END DO
     597             :          END IF
     598             :       END IF
     599             : 
     600          10 :    END SUBROUTINE mao_reference_basis
     601             : 
     602             : ! **************************************************************************************************
     603             : !> \brief Analyze the MAO basis, projection on angular functions
     604             : !> \param mao_coef ...
     605             : !> \param matrix_smm ...
     606             : !> \param mao_basis_set_list ...
     607             : !> \param particle_set ...
     608             : !> \param qs_kind_set ...
     609             : !> \param unit_nr ...
     610             : !> \param para_env ...
     611             : !> \par History
     612             : !>      07.2016 created [JGH]
     613             : ! **************************************************************************************************
     614          10 :    SUBROUTINE mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, &
     615             :                                  qs_kind_set, unit_nr, para_env)
     616             : 
     617             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mao_coef, matrix_smm
     618             :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: mao_basis_set_list
     619             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     620             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     621             :       INTEGER, INTENT(IN)                                :: unit_nr
     622             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     623             : 
     624             :       CHARACTER(len=2)                                   :: element_symbol
     625             :       INTEGER                                            :: ia, iab, iatom, ikind, iset, ishell, &
     626             :                                                             ispin, l, lmax, lshell, m, ma, na, &
     627             :                                                             natom, nspin
     628             :       LOGICAL                                            :: found
     629          10 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cmask, vec1, vec2
     630          10 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: weight
     631          10 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: block, cmao
     632             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
     633             : 
     634             :       ! Analyze the MAO basis
     635          10 :       IF (unit_nr > 0) THEN
     636           5 :          WRITE (unit_nr, "(/,A)") " Analyze angular momentum character of MAOs "
     637             :          WRITE (unit_nr, "(T7,A,T15,A,T20,A,T40,A,T50,A,T60,A,T70,A,T80,A)") &
     638           5 :             "ATOM", "Spin", "MAO", "S", "P", "D", "F", "G"
     639             :       END IF
     640          10 :       lmax = 4 ! analyze up to g-functions
     641          10 :       natom = SIZE(particle_set)
     642          10 :       nspin = SIZE(mao_coef)
     643          58 :       DO iatom = 1, natom
     644             :          CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     645          48 :                               element_symbol=element_symbol, kind_number=ikind)
     646          48 :          basis_set => mao_basis_set_list(ikind)%gto_basis_set
     647          48 :          CALL get_qs_kind(qs_kind_set(ikind), mao=na)
     648          48 :          CALL get_gto_basis_set(basis_set, nsgf=ma)
     649         336 :          ALLOCATE (cmask(ma), vec1(ma), vec2(ma), weight(0:lmax, na))
     650         624 :          weight = 0.0_dp
     651             :          CALL dbcsr_get_block_p(matrix=matrix_smm(1)%matrix, row=iatom, col=iatom, &
     652          48 :                                 block=block, found=found)
     653         102 :          DO ispin = 1, nspin
     654             :             CALL dbcsr_get_block_p(matrix=mao_coef(ispin)%matrix, row=iatom, col=iatom, &
     655          54 :                                    block=cmao, found=found)
     656          54 :             IF (found) THEN
     657         162 :                DO l = 0, lmax
     658        2445 :                   cmask = 0.0_dp
     659         135 :                   iab = 0
     660         675 :                   DO iset = 1, basis_set%nset
     661        1635 :                      DO ishell = 1, basis_set%nshell(iset)
     662         960 :                         lshell = basis_set%l(ishell, iset)
     663        3810 :                         DO m = -lshell, lshell
     664        2310 :                            iab = iab + 1
     665        3270 :                            IF (l == lshell) cmask(iab) = 1.0_dp
     666             :                         END DO
     667             :                      END DO
     668             :                   END DO
     669         432 :                   DO ia = 1, na
     670        6450 :                      vec1(1:ma) = cmask*cmao(1:ma, ia)
     671      383430 :                      vec2(1:ma) = MATMUL(block, vec1)
     672        6585 :                      weight(l, ia) = SUM(vec1(1:ma)*vec2(1:ma))
     673             :                   END DO
     674             :                END DO
     675             :             END IF
     676          54 :             CALL para_env%sum(weight)
     677         156 :             IF (unit_nr > 0) THEN
     678          81 :                DO ia = 1, na
     679          81 :                   IF (ispin == 1 .AND. ia == 1) THEN
     680             :                      WRITE (unit_nr, "(i6,T9,A2,T17,i2,T20,i3,T31,5F10.4)") &
     681          24 :                         iatom, element_symbol, ispin, ia, weight(0:lmax, ia)
     682             :                   ELSE
     683          30 :                      WRITE (unit_nr, "(T17,i2,T20,i3,T31,5F10.4)") ispin, ia, weight(0:lmax, ia)
     684             :                   END IF
     685             :                END DO
     686             :             END IF
     687             :          END DO
     688         154 :          DEALLOCATE (cmask, weight, vec1, vec2)
     689             :       END DO
     690          20 :    END SUBROUTINE mao_basis_analysis
     691             : 
     692             : ! **************************************************************************************************
     693             : !> \brief Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap
     694             : !> \param matrix_q ...
     695             : !> \param matrix_p ...
     696             : !> \param matrix_s ...
     697             : !> \param matrix_smm ...
     698             : !> \param matrix_smo ...
     699             : !> \param smm_list ...
     700             : !> \param electra ...
     701             : !> \param eps_filter ...
     702             : !> \param nimages ...
     703             : !> \param kpoints ...
     704             : !> \param matrix_ks ...
     705             : !> \param sab_orb ...
     706             : !> \par History
     707             : !>      08.2016 created [JGH]
     708             : ! **************************************************************************************************
     709          14 :    SUBROUTINE mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, &
     710             :                           electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb)
     711             : 
     712             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_q
     713             :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p, matrix_s
     714             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_smm, matrix_smo
     715             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     716             :          POINTER                                         :: smm_list
     717             :       REAL(KIND=dp), DIMENSION(2), INTENT(OUT)           :: electra
     718             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
     719             :       INTEGER, INTENT(IN), OPTIONAL                      :: nimages
     720             :       TYPE(kpoint_type), OPTIONAL, POINTER               :: kpoints
     721             :       TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
     722             :          POINTER                                         :: matrix_ks
     723             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     724             :          OPTIONAL, POINTER                               :: sab_orb
     725             : 
     726             :       INTEGER                                            :: im, ispin, nim, nocc, norb, nspin
     727          14 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     728             :       REAL(KIND=dp)                                      :: elex, xkp(3)
     729             :       TYPE(dbcsr_type)                                   :: ksmat, pmat, smat, tmat
     730             : 
     731          14 :       nim = 1
     732          14 :       IF (PRESENT(nimages)) nim = nimages
     733           0 :       IF (nim > 1) THEN
     734           0 :          CPASSERT(PRESENT(kpoints))
     735           0 :          CPASSERT(PRESENT(matrix_ks))
     736           0 :          CPASSERT(PRESENT(sab_orb))
     737             :       END IF
     738             : 
     739             :       ! Reference
     740          14 :       nspin = SIZE(matrix_p, 1)
     741          30 :       DO ispin = 1, nspin
     742          16 :          electra(ispin) = 0.0_dp
     743          46 :          DO im = 1, nim
     744          16 :             CALL dbcsr_dot(matrix_p(ispin, im)%matrix, matrix_s(1, im)%matrix, elex)
     745          32 :             electra(ispin) = electra(ispin) + elex
     746             :          END DO
     747             :       END DO
     748             : 
     749             :       ! Q matrix
     750          14 :       NULLIFY (matrix_q)
     751          14 :       CALL dbcsr_allocate_matrix_set(matrix_q, nspin)
     752          30 :       DO ispin = 1, nspin
     753          16 :          ALLOCATE (matrix_q(ispin)%matrix)
     754          16 :          CALL dbcsr_create(matrix_q(ispin)%matrix, template=matrix_smm(1)%matrix)
     755          30 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_q(ispin)%matrix, smm_list)
     756             :       END DO
     757             :       ! temp matrix
     758          14 :       CALL dbcsr_create(tmat, template=matrix_smo(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
     759             :       ! Q=APA(T)
     760          30 :       DO ispin = 1, nspin
     761          30 :          IF (nim == 1) THEN
     762             :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, matrix_p(ispin, 1)%matrix, &
     763          16 :                                 0.0_dp, tmat, filter_eps=eps_filter)
     764             :             CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
     765          16 :                                 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
     766             :          ELSE
     767             :             ! k-points
     768           0 :             CALL dbcsr_create(pmat, template=matrix_s(1, 1)%matrix)
     769           0 :             CALL dbcsr_create(smat, template=matrix_s(1, 1)%matrix)
     770           0 :             CALL dbcsr_create(ksmat, template=matrix_s(1, 1)%matrix)
     771           0 :             CALL cp_dbcsr_alloc_block_from_nbl(pmat, sab_orb)
     772           0 :             CALL cp_dbcsr_alloc_block_from_nbl(smat, sab_orb)
     773           0 :             CALL cp_dbcsr_alloc_block_from_nbl(ksmat, sab_orb)
     774           0 :             NULLIFY (cell_to_index)
     775           0 :             CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     776             :             ! calculate density matrix at gamma point
     777           0 :             xkp = 0.0_dp
     778             :             ! transform KS and S matrices to the gamma point
     779           0 :             CALL dbcsr_set(ksmat, 0.0_dp)
     780             :             CALL rskp_transform(rmatrix=ksmat, rsmat=matrix_ks, ispin=ispin, &
     781           0 :                                 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
     782           0 :             CALL dbcsr_set(smat, 0.0_dp)
     783             :             CALL rskp_transform(rmatrix=smat, rsmat=matrix_s, ispin=1, &
     784           0 :                                 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb)
     785           0 :             norb = NINT(electra(ispin))
     786           0 :             nocc = MOD(2, nspin) + 1
     787           0 :             CALL calculate_p_gamma(pmat, ksmat, smat, kpoints, norb, REAL(nocc, KIND=dp))
     788             :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, pmat, &
     789           0 :                                 0.0_dp, tmat, filter_eps=eps_filter)
     790             :             CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, &
     791           0 :                                 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter)
     792           0 :             CALL dbcsr_release(pmat)
     793           0 :             CALL dbcsr_release(smat)
     794           0 :             CALL dbcsr_release(ksmat)
     795             :          END IF
     796             :       END DO
     797             :       ! free temp matrix
     798          14 :       CALL dbcsr_release(tmat)
     799             : 
     800          14 :    END SUBROUTINE mao_build_q
     801             : 
     802         864 : END MODULE mao_methods

Generated by: LCOV version 1.15