LCOV - code coverage report
Current view: top level - src - qs_density_matrices.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.8 % 155 150
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 8 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 collects routines that calculate density matrices
      10              : !> \note
      11              : !>      first version : most routines imported
      12              : !> \author JGH (2020-01)
      13              : ! **************************************************************************************************
      14              : MODULE qs_density_matrices
      15              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      16              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      17              :                                               dbcsr_multiply,&
      18              :                                               dbcsr_release,&
      19              :                                               dbcsr_set,&
      20              :                                               dbcsr_type
      21              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_scale_by_vector
      22              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      23              :                                               copy_fm_to_dbcsr,&
      24              :                                               cp_dbcsr_plus_fm_fm_t,&
      25              :                                               cp_dbcsr_sm_fm_multiply
      26              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      27              :                                               cp_fm_scale_and_add,&
      28              :                                               cp_fm_symm,&
      29              :                                               cp_fm_transpose,&
      30              :                                               cp_fm_uplo_to_full
      31              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      32              :                                               cp_fm_struct_release,&
      33              :                                               cp_fm_struct_type
      34              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      35              :                                               cp_fm_get_info,&
      36              :                                               cp_fm_release,&
      37              :                                               cp_fm_to_fm,&
      38              :                                               cp_fm_type
      39              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      40              :                                               cp_logger_get_default_unit_nr,&
      41              :                                               cp_logger_type
      42              :    USE kinds,                           ONLY: dp
      43              :    USE message_passing,                 ONLY: mp_para_env_type
      44              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      45              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      46              :                                               mo_set_type
      47              : #include "./base/base_uses.f90"
      48              : 
      49              :    IMPLICIT NONE
      50              :    PRIVATE
      51              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_density_matrices'
      52              : 
      53              :    PUBLIC :: calculate_density_matrix
      54              :    PUBLIC :: calculate_w_matrix, calculate_w_matrix_ot
      55              :    PUBLIC :: calculate_wz_matrix, calculate_whz_matrix
      56              :    PUBLIC :: calculate_wx_matrix, calculate_xwx_matrix
      57              : 
      58              :    INTERFACE calculate_density_matrix
      59              :       MODULE PROCEDURE calculate_dm_sparse
      60              :    END INTERFACE
      61              : 
      62              :    INTERFACE calculate_w_matrix
      63              :       MODULE PROCEDURE calculate_w_matrix_1, calculate_w_matrix_roks
      64              :    END INTERFACE
      65              : 
      66              : CONTAINS
      67              : 
      68              : ! **************************************************************************************************
      69              : !> \brief   Calculate the density matrix
      70              : !> \param mo_set ...
      71              : !> \param density_matrix ...
      72              : !> \param use_dbcsr ...
      73              : !> \param retain_sparsity ...
      74              : !> \date    06.2002
      75              : !> \par History
      76              : !>       - Fractional occupied orbitals (MK)
      77              : !> \author  Joost VandeVondele
      78              : !> \version 1.0
      79              : ! **************************************************************************************************
      80       205408 :    SUBROUTINE calculate_dm_sparse(mo_set, density_matrix, use_dbcsr, retain_sparsity)
      81              : 
      82              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
      83              :       TYPE(dbcsr_type), POINTER                          :: density_matrix
      84              :       LOGICAL, INTENT(IN), OPTIONAL                      :: use_dbcsr, retain_sparsity
      85              : 
      86              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dm_sparse'
      87              : 
      88              :       INTEGER                                            :: handle
      89              :       LOGICAL                                            :: my_retain_sparsity, my_use_dbcsr
      90              :       REAL(KIND=dp)                                      :: alpha
      91              :       TYPE(cp_fm_type)                                   :: fm_tmp
      92              :       TYPE(dbcsr_type)                                   :: dbcsr_tmp
      93              : 
      94       205408 :       CALL timeset(routineN, handle)
      95              : 
      96       205408 :       my_use_dbcsr = .FALSE.
      97       205408 :       IF (PRESENT(use_dbcsr)) my_use_dbcsr = use_dbcsr
      98       205408 :       my_retain_sparsity = .TRUE.
      99       205408 :       IF (PRESENT(retain_sparsity)) my_retain_sparsity = retain_sparsity
     100       205408 :       IF (my_use_dbcsr) THEN
     101        82479 :          IF (.NOT. ASSOCIATED(mo_set%mo_coeff_b)) THEN
     102            0 :             CPABORT("mo_coeff_b NOT ASSOCIATED")
     103              :          END IF
     104              :       END IF
     105              : 
     106       205408 :       CALL dbcsr_set(density_matrix, 0.0_dp)
     107              : 
     108       205408 :       IF (.NOT. mo_set%uniform_occupation) THEN ! not all orbitals 1..homo are equally occupied
     109        16180 :          IF (my_use_dbcsr) THEN
     110            0 :             CALL dbcsr_copy(dbcsr_tmp, mo_set%mo_coeff_b)
     111              :             CALL dbcsr_scale_by_vector(dbcsr_tmp, mo_set%occupation_numbers(1:mo_set%homo), &
     112            0 :                                        side='right')
     113              :             CALL dbcsr_multiply("N", "T", 1.0_dp, mo_set%mo_coeff_b, dbcsr_tmp, &
     114              :                                 1.0_dp, density_matrix, retain_sparsity=my_retain_sparsity, &
     115            0 :                                 last_k=mo_set%homo)
     116            0 :             CALL dbcsr_release(dbcsr_tmp)
     117              :          ELSE
     118        16180 :             CALL cp_fm_create(fm_tmp, mo_set%mo_coeff%matrix_struct)
     119        16180 :             CALL cp_fm_to_fm(mo_set%mo_coeff, fm_tmp)
     120        16180 :             CALL cp_fm_column_scale(fm_tmp, mo_set%occupation_numbers(1:mo_set%homo))
     121        16180 :             alpha = 1.0_dp
     122              :             CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
     123              :                                        matrix_v=mo_set%mo_coeff, &
     124              :                                        matrix_g=fm_tmp, &
     125              :                                        ncol=mo_set%homo, &
     126        16180 :                                        alpha=alpha)
     127        16180 :             CALL cp_fm_release(fm_tmp)
     128              :          END IF
     129              :       ELSE
     130       189228 :          IF (my_use_dbcsr) THEN
     131              :             CALL dbcsr_multiply("N", "T", mo_set%maxocc, mo_set%mo_coeff_b, mo_set%mo_coeff_b, &
     132              :                                 1.0_dp, density_matrix, retain_sparsity=my_retain_sparsity, &
     133        82479 :                                 last_k=mo_set%homo)
     134              :          ELSE
     135       106749 :             alpha = mo_set%maxocc
     136              :             CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
     137              :                                        matrix_v=mo_set%mo_coeff, &
     138              :                                        ncol=mo_set%homo, &
     139       106749 :                                        alpha=alpha)
     140              :          END IF
     141              :       END IF
     142              : 
     143       205408 :       CALL timestop(handle)
     144              : 
     145       205408 :    END SUBROUTINE calculate_dm_sparse
     146              : 
     147              : ! **************************************************************************************************
     148              : !> \brief Calculate the W matrix from the MO eigenvectors, MO eigenvalues,
     149              : !>       and the MO occupation numbers. Only works if they are eigenstates
     150              : !> \param mo_set type containing the full matrix of the MO and the eigenvalues
     151              : !> \param w_matrix sparse matrix
     152              : !>        error
     153              : !> \par History
     154              : !>         Creation (03.03.03,MK)
     155              : !>         Modification that computes it as a full block, several times (e.g. 20)
     156              : !>               faster at the cost of some additional memory
     157              : !> \author MK
     158              : ! **************************************************************************************************
     159         3616 :    SUBROUTINE calculate_w_matrix_1(mo_set, w_matrix)
     160              : 
     161              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     162              :       TYPE(dbcsr_type), POINTER                          :: w_matrix
     163              : 
     164              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_1'
     165              : 
     166              :       INTEGER                                            :: handle, imo
     167              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eigocc
     168              :       TYPE(cp_fm_type)                                   :: weighted_vectors
     169              : 
     170         3616 :       CALL timeset(routineN, handle)
     171              : 
     172         3616 :       CALL dbcsr_set(w_matrix, 0.0_dp)
     173         3616 :       CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
     174         3616 :       CALL cp_fm_to_fm(mo_set%mo_coeff, weighted_vectors)
     175              : 
     176              :       ! scale every column with the occupation
     177        10798 :       ALLOCATE (eigocc(mo_set%homo))
     178              : 
     179        43480 :       DO imo = 1, mo_set%homo
     180        43480 :          eigocc(imo) = mo_set%eigenvalues(imo)*mo_set%occupation_numbers(imo)
     181              :       END DO
     182         3616 :       CALL cp_fm_column_scale(weighted_vectors, eigocc)
     183         3616 :       DEALLOCATE (eigocc)
     184              : 
     185              :       CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix, &
     186              :                                  matrix_v=mo_set%mo_coeff, &
     187              :                                  matrix_g=weighted_vectors, &
     188         3616 :                                  ncol=mo_set%homo)
     189              : 
     190         3616 :       CALL cp_fm_release(weighted_vectors)
     191              : 
     192         3616 :       CALL timestop(handle)
     193              : 
     194         7232 :    END SUBROUTINE calculate_w_matrix_1
     195              : 
     196              : ! **************************************************************************************************
     197              : !> \brief Calculate the W matrix from the MO coefs, MO derivs
     198              : !>        could overwrite the mo_derivs for increased memory efficiency
     199              : !> \param mo_set type containing the full matrix of the MO coefs
     200              : !>        mo_deriv:
     201              : !> \param mo_deriv ...
     202              : !> \param w_matrix sparse matrix
     203              : !> \param s_matrix sparse matrix for the overlap
     204              : !>        error
     205              : !> \par History
     206              : !>         Creation (JV)
     207              : !> \author MK
     208              : ! **************************************************************************************************
     209         2655 :    SUBROUTINE calculate_w_matrix_ot(mo_set, mo_deriv, w_matrix, s_matrix)
     210              : 
     211              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     212              :       TYPE(dbcsr_type), POINTER                          :: mo_deriv, w_matrix, s_matrix
     213              : 
     214              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_ot'
     215              :       LOGICAL, PARAMETER                                 :: check_gradient = .FALSE., &
     216              :                                                             do_symm = .FALSE.
     217              : 
     218              :       INTEGER                                            :: handle, iounit, ncol_global, nrow_global
     219         2655 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: occupation_numbers, scaling_factor
     220              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     221              :       TYPE(cp_fm_type)                                   :: gradient, h_block, h_block_t, &
     222              :                                                             weighted_vectors
     223              :       TYPE(cp_logger_type), POINTER                      :: logger
     224              : 
     225         2655 :       CALL timeset(routineN, handle)
     226         2655 :       NULLIFY (fm_struct_tmp)
     227              : 
     228              :       CALL cp_fm_get_info(matrix=mo_set%mo_coeff, &
     229              :                           ncol_global=ncol_global, &
     230         2655 :                           nrow_global=nrow_global)
     231              : 
     232         2655 :       CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
     233              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, &
     234              :                                para_env=mo_set%mo_coeff%matrix_struct%para_env, &
     235         2655 :                                context=mo_set%mo_coeff%matrix_struct%context)
     236         2655 :       CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
     237              :       IF (do_symm) CALL cp_fm_create(h_block_t, fm_struct_tmp, name="h block t")
     238         2655 :       CALL cp_fm_struct_release(fm_struct_tmp)
     239              : 
     240         2655 :       CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
     241         7929 :       ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
     242        41210 :       scaling_factor = 2.0_dp*occupation_numbers
     243         2655 :       CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
     244         2655 :       CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
     245         2655 :       DEALLOCATE (scaling_factor)
     246              : 
     247              :       ! the convention seems to require the half here, the factor of two is presumably taken care of
     248              :       ! internally in qs_core_hamiltonian
     249              :       CALL parallel_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 0.5_dp, &
     250         2655 :                          mo_set%mo_coeff, weighted_vectors, 0.0_dp, h_block)
     251              : 
     252              :       IF (do_symm) THEN
     253              :          ! at the minimum things are anyway symmetric, but numerically it might not be the case
     254              :          ! needs some investigation to find out if using this is better
     255              :          CALL cp_fm_transpose(h_block, h_block_t)
     256              :          CALL cp_fm_scale_and_add(0.5_dp, h_block, 0.5_dp, h_block_t)
     257              :       END IF
     258              : 
     259              :       ! this could overwrite the mo_derivs to save the weighted_vectors
     260              :       CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
     261         2655 :                          mo_set%mo_coeff, h_block, 0.0_dp, weighted_vectors)
     262              : 
     263         2655 :       CALL dbcsr_set(w_matrix, 0.0_dp)
     264              :       CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix, &
     265              :                                  matrix_v=mo_set%mo_coeff, &
     266              :                                  matrix_g=weighted_vectors, &
     267         2655 :                                  ncol=mo_set%homo)
     268              : 
     269              :       IF (check_gradient) THEN
     270              :          CALL cp_fm_create(gradient, mo_set%mo_coeff%matrix_struct, "gradient")
     271              :          CALL cp_dbcsr_sm_fm_multiply(s_matrix, weighted_vectors, &
     272              :                                       gradient, ncol_global)
     273              : 
     274              :          ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
     275              :          scaling_factor = 2.0_dp*occupation_numbers
     276              :          CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
     277              :          CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
     278              :          DEALLOCATE (scaling_factor)
     279              : 
     280              :          logger => cp_get_default_logger()
     281              :          IF (logger%para_env%is_source()) THEN
     282              :             iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     283              :             WRITE (iounit, *) " maxabs difference ", &
     284              :                MAXVAL(ABS(weighted_vectors%local_data - 2.0_dp*gradient%local_data))
     285              :          END IF
     286              :          CALL cp_fm_release(gradient)
     287              :       END IF
     288              : 
     289              :       IF (do_symm) CALL cp_fm_release(h_block_t)
     290         2655 :       CALL cp_fm_release(weighted_vectors)
     291         2655 :       CALL cp_fm_release(h_block)
     292              : 
     293         2655 :       CALL timestop(handle)
     294              : 
     295         5310 :    END SUBROUTINE calculate_w_matrix_ot
     296              : 
     297              : ! **************************************************************************************************
     298              : !> \brief Calculate the energy-weighted density matrix W if ROKS is active.
     299              : !>        The W matrix is returned in matrix_w.
     300              : !> \param mo_set ...
     301              : !> \param matrix_ks ...
     302              : !> \param matrix_p ...
     303              : !> \param matrix_w ...
     304              : !> \author 04.05.06,MK
     305              : ! **************************************************************************************************
     306          260 :    SUBROUTINE calculate_w_matrix_roks(mo_set, matrix_ks, matrix_p, matrix_w)
     307              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     308              :       TYPE(dbcsr_type), POINTER                          :: matrix_ks, matrix_p, matrix_w
     309              : 
     310              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_roks'
     311              : 
     312              :       INTEGER                                            :: handle, nao
     313              :       TYPE(cp_blacs_env_type), POINTER                   :: context
     314              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     315              :       TYPE(cp_fm_type)                                   :: ks, p, work
     316              :       TYPE(cp_fm_type), POINTER                          :: c
     317              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     318              : 
     319           52 :       CALL timeset(routineN, handle)
     320              : 
     321           52 :       NULLIFY (context)
     322           52 :       NULLIFY (fm_struct)
     323           52 :       NULLIFY (para_env)
     324              : 
     325           52 :       CALL get_mo_set(mo_set=mo_set, mo_coeff=c)
     326           52 :       CALL cp_fm_get_info(c, context=context, nrow_global=nao, para_env=para_env)
     327              :       CALL cp_fm_struct_create(fm_struct, context=context, nrow_global=nao, &
     328           52 :                                ncol_global=nao, para_env=para_env)
     329           52 :       CALL cp_fm_create(ks, fm_struct, name="Kohn-Sham matrix")
     330           52 :       CALL cp_fm_create(p, fm_struct, name="Density matrix")
     331           52 :       CALL cp_fm_create(work, fm_struct, name="Work matrix")
     332           52 :       CALL cp_fm_struct_release(fm_struct)
     333           52 :       CALL copy_dbcsr_to_fm(matrix_ks, ks)
     334           52 :       CALL copy_dbcsr_to_fm(matrix_p, p)
     335           52 :       CALL cp_fm_uplo_to_full(p, work)
     336           52 :       CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ks, p, 0.0_dp, work)
     337           52 :       CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, p, work, 0.0_dp, ks)
     338           52 :       CALL dbcsr_set(matrix_w, 0.0_dp)
     339           52 :       CALL copy_fm_to_dbcsr(ks, matrix_w, keep_sparsity=.TRUE.)
     340              : 
     341           52 :       CALL cp_fm_release(work)
     342           52 :       CALL cp_fm_release(p)
     343           52 :       CALL cp_fm_release(ks)
     344              : 
     345           52 :       CALL timestop(handle)
     346              : 
     347           52 :    END SUBROUTINE calculate_w_matrix_roks
     348              : 
     349              : ! **************************************************************************************************
     350              : !> \brief Calculate the response W matrix from the MO eigenvectors, MO eigenvalues,
     351              : !>       and the MO occupation numbers. Only works if they are eigenstates
     352              : !> \param mo_set type containing the full matrix of the MO and the eigenvalues
     353              : !> \param psi1 response orbitals
     354              : !> \param ks_matrix Kohn-Sham sparse matrix
     355              : !> \param w_matrix sparse matrix
     356              : !> \par History
     357              : !>               adapted from calculate_w_matrix_1
     358              : !> \author JGH
     359              : ! **************************************************************************************************
     360         4120 :    SUBROUTINE calculate_wz_matrix(mo_set, psi1, ks_matrix, w_matrix)
     361              : 
     362              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     363              :       TYPE(cp_fm_type), INTENT(IN)                       :: psi1
     364              :       TYPE(dbcsr_type), POINTER                          :: ks_matrix, w_matrix
     365              : 
     366              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_wz_matrix'
     367              : 
     368              :       INTEGER                                            :: handle, ncol, nocc, nrow
     369              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     370              :       TYPE(cp_fm_type)                                   :: ksmat, scrv
     371              : 
     372         1030 :       CALL timeset(routineN, handle)
     373              : 
     374              : !     CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
     375              : !     CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
     376              : !     CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     377              : !                              para_env=mo_set%mo_coeff%matrix_struct%para_env, &
     378              : !                              context=mo_set%mo_coeff%matrix_struct%context)
     379              : !     CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
     380              : !     CALL cp_fm_struct_release(fm_struct_tmp)
     381              : !     CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, ncol)
     382              : !     CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
     383              : !     CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
     384              : !     CALL dbcsr_set(w_matrix, 0.0_dp)
     385              : !     CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, &
     386              : !                                ncol=mo_set%homo, symmetry_mode=1)
     387              : !     CALL cp_fm_release(scrv)
     388              : !     CALL cp_fm_release(ksmat)
     389         1030 :       CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
     390         1030 :       nocc = mo_set%homo
     391         1030 :       CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
     392              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nocc, ncol_global=nocc, &
     393              :                                para_env=mo_set%mo_coeff%matrix_struct%para_env, &
     394         1030 :                                context=mo_set%mo_coeff%matrix_struct%context)
     395         1030 :       CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
     396         1030 :       CALL cp_fm_struct_release(fm_struct_tmp)
     397         1030 :       CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc)
     398         1030 :       CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
     399         1030 :       CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
     400         1030 :       CALL dbcsr_set(w_matrix, 0.0_dp)
     401         1030 :       CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1)
     402         1030 :       CALL cp_fm_release(scrv)
     403         1030 :       CALL cp_fm_release(ksmat)
     404              : 
     405         1030 :       CALL timestop(handle)
     406              : 
     407         1030 :    END SUBROUTINE calculate_wz_matrix
     408              : 
     409              : ! **************************************************************************************************
     410              : !> \brief Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues,
     411              : !>       and the MO occupation numbers. Only works if they are eigenstates
     412              : !> \param c0vec ...
     413              : !> \param hzm ...
     414              : !> \param w_matrix sparse matrix
     415              : !> \param focc ...
     416              : !> \param nocc ...
     417              : !> \par History
     418              : !>               adapted from calculate_w_matrix_1
     419              : !> \author JGH
     420              : ! **************************************************************************************************
     421         6040 :    SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
     422              : 
     423              :       TYPE(cp_fm_type), INTENT(IN)                       :: c0vec
     424              :       TYPE(dbcsr_type), POINTER                          :: hzm, w_matrix
     425              :       REAL(KIND=dp), INTENT(IN)                          :: focc
     426              :       INTEGER, INTENT(IN)                                :: nocc
     427              : 
     428              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_whz_matrix'
     429              : 
     430              :       INTEGER                                            :: handle, nao, norb
     431              :       REAL(KIND=dp)                                      :: falpha
     432              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_struct_mat
     433              :       TYPE(cp_fm_type)                                   :: chcmat, hcvec
     434              : 
     435         1510 :       CALL timeset(routineN, handle)
     436              : 
     437         1510 :       falpha = focc
     438              : 
     439         1510 :       CALL cp_fm_create(hcvec, c0vec%matrix_struct, "hcvec")
     440         1510 :       CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb)
     441         1510 :       CPASSERT(nocc <= norb .AND. nocc > 0)
     442         1510 :       norb = nocc
     443              :       CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
     444         1510 :                                ncol_global=norb, para_env=fm_struct%para_env)
     445         1510 :       CALL cp_fm_create(chcmat, fm_struct_mat)
     446         1510 :       CALL cp_fm_struct_release(fm_struct_mat)
     447              : 
     448         1510 :       CALL cp_dbcsr_sm_fm_multiply(hzm, c0vec, hcvec, norb)
     449         1510 :       CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c0vec, hcvec, 0.0_dp, chcmat)
     450         1510 :       CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, c0vec, chcmat, 0.0_dp, hcvec)
     451              : 
     452         1510 :       CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=hcvec, matrix_g=c0vec, ncol=norb, alpha=falpha)
     453              : 
     454         1510 :       CALL cp_fm_release(hcvec)
     455         1510 :       CALL cp_fm_release(chcmat)
     456              : 
     457         1510 :       CALL timestop(handle)
     458              : 
     459         1510 :    END SUBROUTINE calculate_whz_matrix
     460              : 
     461              : ! **************************************************************************************************
     462              : !> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
     463              : !> \param mos_occ ...
     464              : !> \param xvec ...
     465              : !> \param ks_matrix ...
     466              : !> \param w_matrix ...
     467              : !> \par History
     468              : !>               adapted from calculate_wz_matrix
     469              : !> \author JGH
     470              : ! **************************************************************************************************
     471         2656 :    SUBROUTINE calculate_wx_matrix(mos_occ, xvec, ks_matrix, w_matrix)
     472              : 
     473              :       TYPE(cp_fm_type), INTENT(IN)                       :: mos_occ, xvec
     474              :       TYPE(dbcsr_type), POINTER                          :: ks_matrix, w_matrix
     475              : 
     476              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_wx_matrix'
     477              : 
     478              :       INTEGER                                            :: handle, ncol, nrow
     479              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     480              :       TYPE(cp_fm_type)                                   :: ksmat, scrv
     481              : 
     482          664 :       CALL timeset(routineN, handle)
     483              : 
     484          664 :       CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
     485          664 :       CALL cp_fm_create(scrv, mos_occ%matrix_struct, "scr vectors")
     486              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     487              :                                para_env=mos_occ%matrix_struct%para_env, &
     488          664 :                                context=mos_occ%matrix_struct%context)
     489          664 :       CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
     490          664 :       CALL cp_fm_struct_release(fm_struct_tmp)
     491          664 :       CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mos_occ, scrv, ncol)
     492          664 :       CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mos_occ, scrv, 0.0_dp, ksmat)
     493          664 :       CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, xvec, ksmat, 0.0_dp, scrv)
     494          664 :       CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=xvec, ncol=ncol, symmetry_mode=1)
     495          664 :       CALL cp_fm_release(scrv)
     496          664 :       CALL cp_fm_release(ksmat)
     497              : 
     498          664 :       CALL timestop(handle)
     499              : 
     500          664 :    END SUBROUTINE calculate_wx_matrix
     501              : 
     502              : ! **************************************************************************************************
     503              : !> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
     504              : !> \param mos_occ ...
     505              : !> \param xvec ...
     506              : !> \param s_matrix ...
     507              : !> \param ks_matrix ...
     508              : !> \param w_matrix ...
     509              : !> \param eval ...
     510              : !> \par History
     511              : !>               adapted from calculate_wz_matrix
     512              : !> \author JGH
     513              : ! **************************************************************************************************
     514         2656 :    SUBROUTINE calculate_xwx_matrix(mos_occ, xvec, s_matrix, ks_matrix, w_matrix, eval)
     515              : 
     516              :       TYPE(cp_fm_type), INTENT(IN)                       :: mos_occ, xvec
     517              :       TYPE(dbcsr_type), POINTER                          :: s_matrix, ks_matrix, w_matrix
     518              :       REAL(KIND=dp), INTENT(IN)                          :: eval
     519              : 
     520              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_xwx_matrix'
     521              : 
     522              :       INTEGER                                            :: handle, ncol, nrow
     523              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     524              :       TYPE(cp_fm_type)                                   :: scrv, xsxmat
     525              : 
     526          664 :       CALL timeset(routineN, handle)
     527              : 
     528          664 :       CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
     529          664 :       CALL cp_fm_create(scrv, mos_occ%matrix_struct, "scr vectors")
     530              :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     531              :                                para_env=mos_occ%matrix_struct%para_env, &
     532          664 :                                context=mos_occ%matrix_struct%context)
     533          664 :       CALL cp_fm_create(xsxmat, fm_struct_tmp, name="XSX")
     534          664 :       CALL cp_fm_struct_release(fm_struct_tmp)
     535              : 
     536          664 :       CALL cp_dbcsr_sm_fm_multiply(ks_matrix, xvec, scrv, ncol, 1.0_dp, 0.0_dp)
     537          664 :       CALL cp_dbcsr_sm_fm_multiply(s_matrix, xvec, scrv, ncol, eval, -1.0_dp)
     538          664 :       CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, xvec, scrv, 0.0_dp, xsxmat)
     539              : 
     540          664 :       CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mos_occ, xsxmat, 0.0_dp, scrv)
     541          664 :       CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=mos_occ, ncol=ncol, symmetry_mode=1)
     542              : 
     543          664 :       CALL cp_fm_release(scrv)
     544          664 :       CALL cp_fm_release(xsxmat)
     545              : 
     546          664 :       CALL timestop(handle)
     547              : 
     548          664 :    END SUBROUTINE calculate_xwx_matrix
     549              : 
     550              : END MODULE qs_density_matrices
        

Generated by: LCOV version 2.0-1