LCOV - code coverage report
Current view: top level - src - pao_param_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 129 130 99.2 %
Date: 2024-03-28 07:31:50 Functions: 4 4 100.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 Common routines for PAO parametrizations.
      10             : !> \author Ole Schuett
      11             : ! **************************************************************************************************
      12             : MODULE pao_param_methods
      13             :    USE cp_control_types,                ONLY: dft_control_type
      14             :    USE cp_log_handling,                 ONLY: cp_to_string
      15             :    USE dbcsr_api,                       ONLY: &
      16             :         dbcsr_add, dbcsr_complete_redistribute, dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, &
      17             :         dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
      18             :         dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
      19             :         dbcsr_reserve_diag_blocks, dbcsr_scale, dbcsr_type
      20             :    USE dm_ls_scf_qs,                    ONLY: matrix_decluster
      21             :    USE dm_ls_scf_types,                 ONLY: ls_mstruct_type,&
      22             :                                               ls_scf_env_type
      23             :    USE kinds,                           ONLY: dp
      24             :    USE message_passing,                 ONLY: mp_comm_type
      25             :    USE pao_types,                       ONLY: pao_env_type
      26             :    USE qs_environment_types,            ONLY: get_qs_env,&
      27             :                                               qs_environment_type
      28             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      29             :                                               qs_rho_type
      30             : #include "./base/base_uses.f90"
      31             : 
      32             :    IMPLICIT NONE
      33             : 
      34             :    PRIVATE
      35             : 
      36             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_methods'
      37             : 
      38             :    PUBLIC :: pao_calc_grad_lnv_wrt_U, pao_calc_AB_from_U, pao_calc_grad_lnv_wrt_AB
      39             : 
      40             : CONTAINS
      41             : 
      42             : ! **************************************************************************************************
      43             : !> \brief Helper routine, calculates partial derivative dE/dU
      44             : !> \param qs_env ...
      45             : !> \param ls_scf_env ...
      46             : !> \param matrix_M_diag the derivate wrt U, matrix uses pao%diag_distribution
      47             : ! **************************************************************************************************
      48        2432 :    SUBROUTINE pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M_diag)
      49             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      50             :       TYPE(ls_scf_env_type), TARGET                      :: ls_scf_env
      51             :       TYPE(dbcsr_type)                                   :: matrix_M_diag
      52             : 
      53             :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_grad_lnv_wrt_U'
      54             : 
      55             :       INTEGER                                            :: handle
      56             :       REAL(KIND=dp)                                      :: filter_eps
      57        2432 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      58             :       TYPE(dbcsr_type)                                   :: matrix_M, matrix_Ma, matrix_Mb, matrix_NM
      59             :       TYPE(ls_mstruct_type), POINTER                     :: ls_mstruct
      60             :       TYPE(pao_env_type), POINTER                        :: pao
      61             : 
      62        2432 :       CALL timeset(routineN, handle)
      63             : 
      64        2432 :       ls_mstruct => ls_scf_env%ls_mstruct
      65        2432 :       pao => ls_scf_env%pao_env
      66        2432 :       filter_eps = ls_scf_env%eps_filter
      67        2432 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
      68             : 
      69        2432 :       CALL pao_calc_grad_lnv_wrt_AB(qs_env, ls_scf_env, matrix_Ma, matrix_Mb)
      70             : 
      71             :       ! Calculation uses distr. of matrix_s, afterwards we redistribute to pao%diag_distribution.
      72        2432 :       CALL dbcsr_create(matrix_M, template=matrix_s(1)%matrix, matrix_type="N")
      73        2432 :       CALL dbcsr_reserve_diag_blocks(matrix_M)
      74             : 
      75        2432 :       CALL dbcsr_create(matrix_NM, template=ls_mstruct%matrix_A, matrix_type="N")
      76             : 
      77             :       CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_N_inv, matrix_Ma, &
      78        2432 :                           1.0_dp, matrix_NM, filter_eps=filter_eps)
      79             : 
      80             :       CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_N, matrix_Mb, &
      81        2432 :                           1.0_dp, matrix_NM, filter_eps=filter_eps)
      82             : 
      83             :       CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_NM, pao%matrix_Y, &
      84        2432 :                           1.0_dp, matrix_M, filter_eps=filter_eps)
      85             : 
      86             :       !---------------------------------------------------------------------------
      87             :       ! redistribute using pao%diag_distribution
      88             :       CALL dbcsr_create(matrix_M_diag, &
      89             :                         name="PAO matrix_M", &
      90             :                         matrix_type="N", &
      91             :                         dist=pao%diag_distribution, &
      92        2432 :                         template=matrix_s(1)%matrix)
      93        2432 :       CALL dbcsr_reserve_diag_blocks(matrix_M_diag)
      94        2432 :       CALL dbcsr_complete_redistribute(matrix_M, matrix_M_diag)
      95             : 
      96             :       !---------------------------------------------------------------------------
      97             :       ! cleanup:
      98        2432 :       CALL dbcsr_release(matrix_M)
      99        2432 :       CALL dbcsr_release(matrix_Ma)
     100        2432 :       CALL dbcsr_release(matrix_Mb)
     101        2432 :       CALL dbcsr_release(matrix_NM)
     102             : 
     103        2432 :       CALL timestop(handle)
     104        2432 :    END SUBROUTINE pao_calc_grad_lnv_wrt_U
     105             : 
     106             : ! **************************************************************************************************
     107             : !> \brief Takes current matrix_X and calculates the matrices A and B.
     108             : !> \param pao ...
     109             : !> \param qs_env ...
     110             : !> \param ls_scf_env ...
     111             : !> \param matrix_U_diag ...
     112             : ! **************************************************************************************************
     113       13070 :    SUBROUTINE pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U_diag)
     114             :       TYPE(pao_env_type), POINTER                        :: pao
     115             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     116             :       TYPE(ls_scf_env_type), TARGET                      :: ls_scf_env
     117             :       TYPE(dbcsr_type)                                   :: matrix_U_diag
     118             : 
     119             :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_AB_from_U'
     120             : 
     121             :       INTEGER                                            :: acol, arow, handle, iatom
     122             :       LOGICAL                                            :: found
     123       13070 :       REAL(dp), DIMENSION(:, :), POINTER                 :: block_A, block_B, block_N, block_N_inv, &
     124       13070 :                                                             block_U, block_Y
     125             :       TYPE(dbcsr_iterator_type)                          :: iter
     126       13070 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     127             :       TYPE(dbcsr_type)                                   :: matrix_U
     128             :       TYPE(ls_mstruct_type), POINTER                     :: ls_mstruct
     129             : 
     130       13070 :       CALL timeset(routineN, handle)
     131       13070 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     132       13070 :       ls_mstruct => ls_scf_env%ls_mstruct
     133             : 
     134             :       ! --------------------------------------------------------------------------------------------
     135             :       ! sanity check matrix U
     136       13070 :       CALL pao_assert_unitary(pao, matrix_U_diag)
     137             : 
     138             :       ! --------------------------------------------------------------------------------------------
     139             :       ! redistribute matrix_U_diag from diag_distribution to distribution of matrix_s
     140       13070 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     141       13070 :       CALL dbcsr_create(matrix_U, matrix_type="N", template=matrix_s(1)%matrix)
     142       13070 :       CALL dbcsr_reserve_diag_blocks(matrix_U)
     143       13070 :       CALL dbcsr_complete_redistribute(matrix_U_diag, matrix_U)
     144             : 
     145             :       ! --------------------------------------------------------------------------------------------
     146             :       ! calculate matrix A and B from matrix U
     147             :       ! Multiplying diagonal matrices is a local operation.
     148             :       ! To take advantage of this we're using an iterator instead of calling dbcsr_multiply().
     149             : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,ls_mstruct,matrix_U) &
     150       13070 : !$OMP PRIVATE(iter,arow,acol,iatom,block_U,block_Y,block_A,block_B,block_N,block_N_inv,found)
     151             :       CALL dbcsr_iterator_start(iter, matrix_U)
     152             :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     153             :          CALL dbcsr_iterator_next_block(iter, arow, acol, block_U)
     154             :          iatom = arow; CPASSERT(arow == acol)
     155             : 
     156             :          CALL dbcsr_get_block_p(matrix=pao%matrix_Y, row=iatom, col=iatom, block=block_Y, found=found)
     157             :          CPASSERT(ASSOCIATED(block_Y))
     158             : 
     159             :          CALL dbcsr_get_block_p(matrix=ls_mstruct%matrix_A, row=iatom, col=iatom, block=block_A, found=found)
     160             :          CALL dbcsr_get_block_p(matrix=pao%matrix_N_inv, row=iatom, col=iatom, block=block_N_inv, found=found)
     161             :          CPASSERT(ASSOCIATED(block_A) .AND. ASSOCIATED(block_N_inv))
     162             : 
     163             :          CALL dbcsr_get_block_p(matrix=ls_mstruct%matrix_B, row=iatom, col=iatom, block=block_B, found=found)
     164             :          CALL dbcsr_get_block_p(matrix=pao%matrix_N, row=iatom, col=iatom, block=block_N, found=found)
     165             :          CPASSERT(ASSOCIATED(block_B) .AND. ASSOCIATED(block_N))
     166             : 
     167             :          block_A = MATMUL(MATMUL(block_N_inv, block_U), block_Y)
     168             :          block_B = MATMUL(MATMUL(block_N, block_U), block_Y)
     169             : 
     170             :       END DO
     171             :       CALL dbcsr_iterator_stop(iter)
     172             : !$OMP END PARALLEL
     173             : 
     174       13070 :       CALL dbcsr_release(matrix_U)
     175             : 
     176       13070 :       CALL timestop(handle)
     177       13070 :    END SUBROUTINE pao_calc_AB_from_U
     178             : 
     179             : ! **************************************************************************************************
     180             : !> \brief Debugging routine, check unitaryness of U
     181             : !> \param pao ...
     182             : !> \param matrix_U ...
     183             : ! **************************************************************************************************
     184       18018 :    SUBROUTINE pao_assert_unitary(pao, matrix_U)
     185             :       TYPE(pao_env_type), POINTER                        :: pao
     186             :       TYPE(dbcsr_type)                                   :: matrix_U
     187             : 
     188             :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_assert_unitary'
     189             : 
     190             :       INTEGER                                            :: acol, arow, group_handle, handle, i, &
     191             :                                                             iatom, M, N
     192       13070 :       INTEGER, DIMENSION(:), POINTER                     :: blk_sizes_pao, blk_sizes_pri
     193             :       REAL(dp)                                           :: delta_max
     194       13070 :       REAL(dp), DIMENSION(:, :), POINTER                 :: block_test, tmp1, tmp2
     195             :       TYPE(dbcsr_iterator_type)                          :: iter
     196             :       TYPE(mp_comm_type)                                 :: group
     197             : 
     198       10596 :       IF (pao%check_unitary_tol < 0.0_dp) RETURN ! no checking
     199             : 
     200        2474 :       CALL timeset(routineN, handle)
     201        2474 :       delta_max = 0.0_dp
     202             : 
     203        2474 :       CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
     204             : 
     205             : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_U,blk_sizes_pri,blk_sizes_pao,delta_max) &
     206        2474 : !$OMP PRIVATE(iter,arow,acol,iatom,N,M,block_test,tmp1,tmp2)
     207             :       CALL dbcsr_iterator_start(iter, matrix_U)
     208             :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     209             :          CALL dbcsr_iterator_next_block(iter, arow, acol, block_test)
     210             :          iatom = arow; CPASSERT(arow == acol)
     211             :          N = blk_sizes_pri(iatom) ! size of primary basis
     212             :          M = blk_sizes_pao(iatom) ! size of pao basis
     213             : 
     214             :          ! we only need the upper left "PAO-corner" to be unitary
     215             :          ALLOCATE (tmp1(N, M), tmp2(M, M))
     216             :          tmp1 = block_test(:, 1:M)
     217             :          tmp2 = MATMUL(TRANSPOSE(tmp1), tmp1)
     218             :          DO i = 1, M
     219             :             tmp2(i, i) = tmp2(i, i) - 1.0_dp
     220             :          END DO
     221             : 
     222             : !$OMP ATOMIC
     223             :          delta_max = MAX(delta_max, MAXVAL(ABS(tmp2)))
     224             : 
     225             :          DEALLOCATE (tmp1, tmp2)
     226             :       END DO
     227             :       CALL dbcsr_iterator_stop(iter)
     228             : !$OMP END PARALLEL
     229             : 
     230        2474 :       CALL dbcsr_get_info(matrix_U, group=group_handle)
     231        2474 :       CALL group%set_handle(group_handle)
     232             : 
     233        2474 :       CALL group%max(delta_max)
     234        2474 :       IF (pao%iw > 0) WRITE (pao%iw, *) 'PAO| checked unitaryness, max delta:', delta_max
     235        2474 :       IF (delta_max > pao%check_unitary_tol) &
     236           0 :          CPABORT("Found bad unitaryness:"//cp_to_string(delta_max))
     237             : 
     238        2474 :       CALL timestop(handle)
     239       13070 :    END SUBROUTINE pao_assert_unitary
     240             : 
     241             : ! **************************************************************************************************
     242             : !> \brief Helper routine, calculates partial derivative dE/dA and dE/dB.
     243             : !>        As energy functional serves the definition by LNV (Li, Nunes, Vanderbilt).
     244             : !> \param qs_env ...
     245             : !> \param ls_scf_env ...
     246             : !> \param matrix_Ma the derivate wrt A, matrix uses s_matrix-distribution.
     247             : !> \param matrix_Mb the derivate wrt B, matrix uses s_matrix-distribution.
     248             : ! **************************************************************************************************
     249        2664 :    SUBROUTINE pao_calc_grad_lnv_wrt_AB(qs_env, ls_scf_env, matrix_Ma, matrix_Mb)
     250             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     251             :       TYPE(ls_scf_env_type), TARGET                      :: ls_scf_env
     252             :       TYPE(dbcsr_type)                                   :: matrix_Ma, matrix_Mb
     253             : 
     254             :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_grad_lnv_wrt_AB'
     255             : 
     256             :       INTEGER                                            :: handle, nspin
     257        2664 :       INTEGER, DIMENSION(:), POINTER                     :: pao_blk_sizes
     258             :       REAL(KIND=dp)                                      :: filter_eps
     259        2664 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s, rho_ao
     260             :       TYPE(dbcsr_type) :: matrix_HB, matrix_HPS, matrix_M, matrix_M1, matrix_M1_dc, matrix_M2, &
     261             :          matrix_M2_dc, matrix_M3, matrix_M3_dc, matrix_PA, matrix_PH, matrix_PHP, matrix_PSP, &
     262             :          matrix_SB, matrix_SP
     263             :       TYPE(dft_control_type), POINTER                    :: dft_control
     264             :       TYPE(ls_mstruct_type), POINTER                     :: ls_mstruct
     265             :       TYPE(pao_env_type), POINTER                        :: pao
     266             :       TYPE(qs_rho_type), POINTER                         :: rho
     267             : 
     268        2664 :       CALL timeset(routineN, handle)
     269             : 
     270        2664 :       ls_mstruct => ls_scf_env%ls_mstruct
     271        2664 :       pao => ls_scf_env%pao_env
     272             : 
     273             :       CALL get_qs_env(qs_env, &
     274             :                       rho=rho, &
     275             :                       matrix_ks=matrix_ks, &
     276             :                       matrix_s=matrix_s, &
     277        2664 :                       dft_control=dft_control)
     278        2664 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     279        2664 :       nspin = dft_control%nspins
     280        2664 :       filter_eps = ls_scf_env%eps_filter
     281             : 
     282        2664 :       CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
     283             : 
     284        2664 :       IF (nspin /= 1) CPABORT("open shell not yet implemented")
     285             :       !TODO: handle openshell case properly
     286             : 
     287             :       ! Notation according to equation (4.6) on page 50 from:
     288             :       ! https://dx.doi.org/10.3929%2Fethz-a-010819495
     289             : 
     290             :       !---------------------------------------------------------------------------
     291             :       ! calculate need products in pao basis
     292        2664 :       CALL dbcsr_create(matrix_PH, template=ls_scf_env%matrix_s, matrix_type="N")
     293             :       CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_p(1), ls_scf_env%matrix_ks(1), &
     294        2664 :                           0.0_dp, matrix_PH, filter_eps=filter_eps)
     295             : 
     296        2664 :       CALL dbcsr_create(matrix_PHP, template=ls_scf_env%matrix_s, matrix_type="N")
     297             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_PH, ls_scf_env%matrix_p(1), &
     298        2664 :                           0.0_dp, matrix_PHP, filter_eps=filter_eps)
     299             : 
     300        2664 :       CALL dbcsr_create(matrix_SP, template=ls_scf_env%matrix_s, matrix_type="N")
     301             :       CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s, ls_scf_env%matrix_p(1), &
     302        2664 :                           0.0_dp, matrix_SP, filter_eps=filter_eps)
     303             : 
     304        2664 :       IF (nspin == 1) CALL dbcsr_scale(matrix_SP, 0.5_dp)
     305             : 
     306        2664 :       CALL dbcsr_create(matrix_HPS, template=ls_scf_env%matrix_s, matrix_type="N")
     307             :       CALL dbcsr_multiply("N", "T", 1.0_dp, ls_scf_env%matrix_ks(1), matrix_SP, &
     308        2664 :                           0.0_dp, matrix_HPS, filter_eps=filter_eps)
     309             : 
     310        2664 :       CALL dbcsr_create(matrix_PSP, template=ls_scf_env%matrix_s, matrix_type="N")
     311             :       CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_p(1), matrix_SP, &
     312        2664 :                           0.0_dp, matrix_PSP, filter_eps=filter_eps)
     313             : 
     314             :       !---------------------------------------------------------------------------
     315             :       ! M1 = dE_lnv / dP_pao
     316        2664 :       CALL dbcsr_create(matrix_M1, template=ls_scf_env%matrix_s, matrix_type="N")
     317             : 
     318             :       CALL dbcsr_multiply("N", "T", 3.0_dp, ls_scf_env%matrix_ks(1), matrix_SP, &
     319        2664 :                           1.0_dp, matrix_M1, filter_eps=filter_eps)
     320             : 
     321             :       CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_SP, ls_scf_env%matrix_ks(1), &
     322        2664 :                           1.0_dp, matrix_M1, filter_eps=filter_eps)
     323             : 
     324             :       CALL dbcsr_multiply("N", "T", -2.0_dp, matrix_HPS, matrix_SP, &
     325        2664 :                           1.0_dp, matrix_M1, filter_eps=filter_eps)
     326             : 
     327             :       CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_SP, matrix_HPS, &
     328        2664 :                           1.0_dp, matrix_M1, filter_eps=filter_eps)
     329             : 
     330             :       CALL dbcsr_multiply("N", "T", -2.0_dp, matrix_SP, matrix_HPS, &
     331        2664 :                           1.0_dp, matrix_M1, filter_eps=filter_eps)
     332             : 
     333             :       ! reverse possible molecular clustering
     334             :       CALL dbcsr_create(matrix_M1_dc, &
     335             :                         template=matrix_s(1)%matrix, &
     336             :                         row_blk_size=pao_blk_sizes, &
     337        2664 :                         col_blk_size=pao_blk_sizes)
     338        2664 :       CALL matrix_decluster(matrix_M1_dc, matrix_M1, ls_mstruct)
     339             : 
     340             :       !---------------------------------------------------------------------------
     341             :       ! M2 = dE_lnv / dH
     342        2664 :       CALL dbcsr_create(matrix_M2, template=ls_scf_env%matrix_s, matrix_type="N")
     343             : 
     344        2664 :       CALL dbcsr_add(matrix_M2, matrix_PSP, 1.0_dp, 3.0_dp)
     345             : 
     346             :       CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_PSP, matrix_SP, &
     347        2664 :                           1.0_dp, matrix_M2, filter_eps=filter_eps)
     348             : 
     349             :       ! reverse possible molecular clustering
     350             :       CALL dbcsr_create(matrix_M2_dc, &
     351             :                         template=matrix_s(1)%matrix, &
     352             :                         row_blk_size=pao_blk_sizes, &
     353        2664 :                         col_blk_size=pao_blk_sizes)
     354        2664 :       CALL matrix_decluster(matrix_M2_dc, matrix_M2, ls_mstruct)
     355             : 
     356             :       !---------------------------------------------------------------------------
     357             :       ! M3 = dE_lnv / dS
     358        2664 :       CALL dbcsr_create(matrix_M3, template=ls_scf_env%matrix_s, matrix_type="N")
     359             : 
     360        2664 :       CALL dbcsr_add(matrix_M3, matrix_PHP, 1.0_dp, 3.0_dp)
     361             : 
     362             :       CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_PHP, matrix_SP, &
     363        2664 :                           1.0_dp, matrix_M3, filter_eps=filter_eps)
     364             : 
     365             :       CALL dbcsr_multiply("N", "T", -2.0_dp, matrix_PSP, matrix_PH, &
     366        2664 :                           1.0_dp, matrix_M3, filter_eps=filter_eps)
     367             : 
     368             :       ! reverse possible molecular clustering
     369             :       CALL dbcsr_create(matrix_M3_dc, &
     370             :                         template=matrix_s(1)%matrix, &
     371             :                         row_blk_size=pao_blk_sizes, &
     372        2664 :                         col_blk_size=pao_blk_sizes)
     373        2664 :       CALL matrix_decluster(matrix_M3_dc, matrix_M3, ls_mstruct)
     374             : 
     375             :       !---------------------------------------------------------------------------
     376             :       ! assemble Ma and Mb
     377             :       ! matrix_Ma = dE_lnv / dA = P * A * M1
     378             :       ! matrix_Mb = dE_lnv / dB = H * B * M2  +  S * B * M3
     379        2664 :       CALL dbcsr_create(matrix_Ma, template=ls_mstruct%matrix_A, matrix_type="N")
     380        2664 :       CALL dbcsr_reserve_diag_blocks(matrix_Ma)
     381        2664 :       CALL dbcsr_create(matrix_Mb, template=ls_mstruct%matrix_B, matrix_type="N")
     382        2664 :       CALL dbcsr_reserve_diag_blocks(matrix_Mb)
     383             : 
     384             :       !---------------------------------------------------------------------------
     385             :       ! combine M1 with matrices from primary basis
     386        2664 :       CALL dbcsr_create(matrix_PA, template=ls_mstruct%matrix_A, matrix_type="N")
     387             :       CALL dbcsr_multiply("N", "N", 1.0_dp, rho_ao(1)%matrix, ls_mstruct%matrix_A, &
     388        2664 :                           0.0_dp, matrix_PA, filter_eps=filter_eps)
     389             : 
     390             :       ! matrix_Ma = P * A * M1
     391             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_PA, matrix_M1_dc, &
     392        2664 :                           0.0_dp, matrix_Ma, filter_eps=filter_eps)
     393             : 
     394             :       !---------------------------------------------------------------------------
     395             :       ! combine M2 with matrices from primary basis
     396        2664 :       CALL dbcsr_create(matrix_HB, template=ls_mstruct%matrix_B, matrix_type="N")
     397             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ks(1)%matrix, ls_mstruct%matrix_B, &
     398        2664 :                           0.0_dp, matrix_HB, filter_eps=filter_eps)
     399             : 
     400             :       ! matrix_Mb = H * B * M2
     401             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_HB, matrix_M2_dc, &
     402        2664 :                           0.0_dp, matrix_Mb, filter_eps=filter_eps)
     403             : 
     404             :       !---------------------------------------------------------------------------
     405             :       ! combine M3 with matrices from primary basis
     406        2664 :       CALL dbcsr_create(matrix_SB, template=ls_mstruct%matrix_B, matrix_type="N")
     407             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s(1)%matrix, ls_mstruct%matrix_B, &
     408        2664 :                           0.0_dp, matrix_SB, filter_eps=filter_eps)
     409             : 
     410        2664 :       IF (nspin == 1) CALL dbcsr_scale(matrix_SB, 0.5_dp)
     411             : 
     412             :       ! matrix_Mb += S * B * M3
     413             :       CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_SB, matrix_M3_dc, &
     414        2664 :                           1.0_dp, matrix_Mb, filter_eps=filter_eps)
     415             : 
     416        2664 :       IF (nspin == 1) CALL dbcsr_scale(matrix_Ma, 2.0_dp)
     417        2664 :       IF (nspin == 1) CALL dbcsr_scale(matrix_Mb, 2.0_dp)
     418             : 
     419             :       !---------------------------------------------------------------------------
     420             :       ! cleanup: TODO release matrices as early as possible
     421        2664 :       CALL dbcsr_release(matrix_PH)
     422        2664 :       CALL dbcsr_release(matrix_PHP)
     423        2664 :       CALL dbcsr_release(matrix_SP)
     424        2664 :       CALL dbcsr_release(matrix_HPS)
     425        2664 :       CALL dbcsr_release(matrix_PSP)
     426        2664 :       CALL dbcsr_release(matrix_M)
     427        2664 :       CALL dbcsr_release(matrix_M1)
     428        2664 :       CALL dbcsr_release(matrix_M2)
     429        2664 :       CALL dbcsr_release(matrix_M3)
     430        2664 :       CALL dbcsr_release(matrix_M1_dc)
     431        2664 :       CALL dbcsr_release(matrix_M2_dc)
     432        2664 :       CALL dbcsr_release(matrix_M3_dc)
     433        2664 :       CALL dbcsr_release(matrix_PA)
     434        2664 :       CALL dbcsr_release(matrix_HB)
     435        2664 :       CALL dbcsr_release(matrix_SB)
     436             : 
     437        2664 :       CALL timestop(handle)
     438        2664 :    END SUBROUTINE pao_calc_grad_lnv_wrt_AB
     439             : 
     440             : END MODULE pao_param_methods

Generated by: LCOV version 1.15