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             : MODULE qs_tddfpt2_smearing_methods
       9             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      10             :    USE cp_control_types,                ONLY: dft_control_type,&
      11             :                                               smeared_type,&
      12             :                                               tddfpt2_control_type
      13             :    USE cp_dbcsr_api,                    ONLY: dbcsr_type
      14             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      15             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      16             :                                               cp_fm_scale_and_add
      17             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      18             :                                               cp_fm_struct_release,&
      19             :                                               cp_fm_struct_type
      20             :    USE cp_fm_types,                     ONLY: cp_fm_copy_general,&
      21             :                                               cp_fm_create,&
      22             :                                               cp_fm_get_info,&
      23             :                                               cp_fm_get_submatrix,&
      24             :                                               cp_fm_release,&
      25             :                                               cp_fm_set_submatrix,&
      26             :                                               cp_fm_to_fm,&
      27             :                                               cp_fm_type
      28             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      29             :                                               cp_logger_get_default_io_unit,&
      30             :                                               cp_logger_type
      31             :    USE fermi_utils,                     ONLY: Fermi
      32             :    USE input_constants,                 ONLY: smear_fermi_dirac
      33             :    USE kinds,                           ONLY: dp
      34             :    USE message_passing,                 ONLY: mp_para_env_type
      35             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      36             :    USE qs_environment_types,            ONLY: get_qs_env,&
      37             :                                               qs_environment_type
      38             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      39             :                                               mo_set_type
      40             :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos
      41             :    USE scf_control_types,               ONLY: scf_control_type,&
      42             :                                               smear_type
      43             : #include "./base/base_uses.f90"
      44             : 
      45             :    IMPLICIT NONE
      46             : 
      47             :    PRIVATE
      48             : 
      49             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_smearing_methods'
      50             : 
      51             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      52             : 
      53             :    PUBLIC :: tddfpt_smeared_occupation, &
      54             :              add_smearing_aterm, compute_fermib, &
      55             :              orthogonalize_smeared_occupation, deallocate_fermi_params
      56             : 
      57             : CONTAINS
      58             : 
      59             : ! **************************************************************************************************
      60             : !> \brief ...
      61             : !> \param qs_env ...
      62             : !> \param gs_mos ...
      63             : !> \param log_unit ...
      64             : ! **************************************************************************************************
      65           2 :    SUBROUTINE tddfpt_smeared_occupation(qs_env, gs_mos, log_unit)
      66             : 
      67             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      68             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
      69             :          INTENT(IN), POINTER                             :: gs_mos
      70             :       INTEGER, INTENT(in)                                :: log_unit
      71             : 
      72             :       CHARACTER(len=*), PARAMETER :: routineN = 'tddfpt_smeared_occupation'
      73             : 
      74             :       INTEGER                                            :: handle, iocc, ispin, nspins
      75           2 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nocc, nvirt
      76           2 :       REAL(kind=dp), DIMENSION(:), POINTER               :: mo_evals, occup
      77           2 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
      78             :       TYPE(scf_control_type), POINTER                    :: scf_control
      79             :       TYPE(smear_type), POINTER                          :: smear
      80             : 
      81           2 :       CALL timeset(routineN, handle)
      82             : 
      83           2 :       nspins = SIZE(gs_mos)
      84             : 
      85           2 :       NULLIFY (mos, scf_control)
      86           2 :       CALL get_qs_env(qs_env, mos=mos, scf_control=scf_control)
      87           2 :       NULLIFY (smear)
      88           2 :       IF (ASSOCIATED(qs_env%scf_control%smear)) THEN
      89           2 :          smear => qs_env%scf_control%smear
      90             :       ELSE
      91           0 :          CPABORT("Smeared input section no longer associated.")
      92             :       END IF
      93             : 
      94             :       IF (debug_this_module) THEN
      95             :          NULLIFY (mo_evals, occup)
      96             :          ALLOCATE (nocc(nspins), nvirt(nspins))
      97             :          DO ispin = 1, nspins
      98             :             CALL get_mo_set(mos(ispin), eigenvalues=mo_evals, occupation_numbers=occup)
      99             :             CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nocc(ispin))
     100             :             CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, ncol_global=nvirt(ispin))
     101             :             IF (log_unit > 0) THEN
     102             :                DO iocc = 1, nocc(ispin)
     103             :                   WRITE (log_unit, '(A,F14.5)') "Occupation numbers", occup(iocc)
     104             :                END DO
     105             :             END IF
     106             :          END DO
     107             :       END IF
     108             : 
     109           2 :       CALL allocate_fermi_params(qs_env, gs_mos)
     110           2 :       CALL compute_fermia(qs_env, gs_mos, log_unit)
     111             : 
     112           2 :       CALL timestop(handle)
     113             : 
     114           2 :    END SUBROUTINE tddfpt_smeared_occupation
     115             : ! **************************************************************************************************
     116             : !> \brief ...
     117             : !> \param qs_env ...
     118             : !> \param gs_mos ...
     119             : !> \param log_unit ...
     120             : ! **************************************************************************************************
     121           2 :    SUBROUTINE compute_fermia(qs_env, gs_mos, log_unit)
     122             : 
     123             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     124             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     125             :          INTENT(IN), POINTER                             :: gs_mos
     126             :       INTEGER, INTENT(IN)                                :: log_unit
     127             : 
     128             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'compute_fermia'
     129             : 
     130             :       INTEGER                                            :: handle, iocc, ispin, nspins
     131           2 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nocc
     132             :       REAL(kind=dp)                                      :: maxvalue, mu
     133           2 :       REAL(kind=dp), DIMENSION(:), POINTER               :: fermia, mo_evals, occup
     134             :       TYPE(dft_control_type), POINTER                    :: dft_control
     135           2 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     136             :       TYPE(scf_control_type), POINTER                    :: scf_control
     137             :       TYPE(smear_type), POINTER                          :: smear
     138             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     139             : 
     140           2 :       CALL timeset(routineN, handle)
     141             : 
     142           2 :       NULLIFY (mos, dft_control, tddfpt_control, scf_control)
     143           2 :       CALL get_qs_env(qs_env, dft_control=dft_control, scf_control=scf_control)
     144           2 :       tddfpt_control => dft_control%tddfpt2_control
     145             : 
     146           2 :       CALL get_qs_env(qs_env, mos=mos)
     147           2 :       nspins = SIZE(mos)
     148             : 
     149           2 :       NULLIFY (smear)
     150           2 :       IF (ASSOCIATED(qs_env%scf_control%smear)) THEN
     151             :          smear => qs_env%scf_control%smear
     152             :       ELSE
     153           0 :          CPABORT("Smeared input section no longer associated.")
     154             :       END IF
     155             : 
     156             :       IF (debug_this_module .AND. (log_unit > 0)) THEN
     157             :          WRITE (log_unit, '(A,F14.5)') "Smearing temperature", smear%electronic_temperature
     158             :       END IF
     159             : 
     160           2 :       NULLIFY (mo_evals, occup)
     161           6 :       ALLOCATE (nocc(nspins))
     162           4 :       DO ispin = 1, nspins
     163           2 :          CALL get_mo_set(mos(ispin), eigenvalues=mo_evals, occupation_numbers=occup, mu=mu)
     164           4 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nocc(ispin))
     165             :       END DO
     166             : 
     167           4 :       DO ispin = 1, nspins
     168           2 :          fermia => tddfpt_control%smeared_occup(ispin)%fermia
     169          12 :          DO iocc = 1, nocc(ispin)
     170           8 :             maxvalue = mu + 3.0_dp*smear%electronic_temperature - mo_evals(iocc)
     171             : 
     172           8 :             fermia(iocc) = MAX(0.0_dp, maxvalue)
     173           2 :             IF (debug_this_module) THEN
     174             :                IF (log_unit > 0) WRITE (log_unit, '(A,F14.5)') "Fermi smearing parameter alpha", fermia(iocc)
     175             :             END IF
     176             :          END DO
     177             :       END DO
     178             : 
     179           2 :       CALL timestop(handle)
     180           4 :    END SUBROUTINE compute_fermia
     181             : ! **************************************************************************************************
     182             : 
     183             : ! **************************************************************************************************
     184             : !> \brief ...
     185             : !> \param qs_env ...
     186             : !> \param Aop_evects ...
     187             : !> \param evects ...
     188             : !> \param S_evects ...
     189             : !> \param mos_occ ...
     190             : !> \param fermia ...
     191             : !> \param matrix_s ...
     192             : ! **************************************************************************************************
     193          72 :    SUBROUTINE add_smearing_aterm(qs_env, Aop_evects, evects, S_evects, mos_occ, fermia, matrix_s)
     194             : 
     195             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     196             :       TYPE(cp_fm_type), INTENT(in)                       :: Aop_evects, evects, S_evects, mos_occ
     197             :       REAL(kind=dp), DIMENSION(:), POINTER               :: fermia
     198             :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     199             : 
     200             :       CHARACTER(len=*), PARAMETER :: routineN = 'add_smearing_aterm'
     201             : 
     202             :       INTEGER                                            :: handle, iounit, nactive, nao
     203             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct, matrix_struct_tmp
     204             :       TYPE(cp_fm_type)                                   :: CCSXvec, CSXvec, Cvec, SCCSXvec
     205             :       TYPE(cp_logger_type), POINTER                      :: logger
     206             :       TYPE(dft_control_type), POINTER                    :: dft_control
     207             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     208             : 
     209          12 :       CALL timeset(routineN, handle)
     210             : 
     211          12 :       NULLIFY (logger)
     212          12 :       logger => cp_get_default_logger()
     213          12 :       iounit = cp_logger_get_default_io_unit(logger)
     214             : 
     215          12 :       NULLIFY (dft_control, tddfpt_control)
     216          12 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     217          12 :       tddfpt_control => dft_control%tddfpt2_control
     218             : 
     219             :       CALL cp_fm_get_info(matrix=evects, matrix_struct=matrix_struct, &
     220          12 :                           nrow_global=nao, ncol_global=nactive)
     221          12 :       CALL cp_fm_create(CCSXvec, matrix_struct)
     222          12 :       CALL cp_fm_create(SCCSXvec, matrix_struct)
     223          12 :       CALL cp_fm_create(Cvec, matrix_struct)
     224             : 
     225          12 :       NULLIFY (matrix_struct_tmp)
     226             :       CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nactive, &
     227          12 :                                ncol_global=nactive, template_fmstruct=matrix_struct)
     228          12 :       CALL cp_fm_create(CSXvec, matrix_struct_tmp)
     229          12 :       CALL cp_fm_struct_release(fmstruct=matrix_struct_tmp)
     230             : 
     231             :       CALL parallel_gemm('T', 'N', nactive, nactive, nao, 1.0_dp, &
     232          12 :                          mos_occ, S_evects, 0.0_dp, CSXvec)
     233          12 :       CALL cp_fm_to_fm(mos_occ, Cvec)
     234             : 
     235          12 :       CALL cp_fm_column_scale(Cvec, fermia)
     236             :       CALL parallel_gemm('N', 'N', nao, nactive, nactive, 1.0_dp, &
     237          12 :                          Cvec, CSXvec, 0.0_dp, CCSXvec)
     238             : 
     239             :       ! alpha S C C^T S X
     240             :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, CCSXvec, &
     241          12 :                                    SCCSXvec, ncol=nactive, alpha=1.0_dp, beta=0.0_dp)
     242          12 :       CALL cp_fm_scale_and_add(1.0_dp, Aop_evects, 1.0_dp, SCCSXvec)
     243             : 
     244          12 :       CALL cp_fm_release(SCCSXvec)
     245          12 :       CALL cp_fm_release(CCSXvec)
     246          12 :       CALL cp_fm_release(CSXvec)
     247          12 :       CALL cp_fm_release(Cvec)
     248             : 
     249          12 :       CALL timestop(handle)
     250          12 :    END SUBROUTINE add_smearing_aterm
     251             : ! **************************************************************************************************
     252             : !> \brief ...
     253             : !> \param qs_env ...
     254             : !> \param gs_mos ...
     255             : !> \param evals ...
     256             : ! **************************************************************************************************
     257          14 :    SUBROUTINE compute_fermib(qs_env, gs_mos, evals)
     258             : 
     259             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     260             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     261             :          INTENT(IN)                                      :: gs_mos
     262             :       REAL(kind=dp), INTENT(in)                          :: evals
     263             : 
     264             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'compute_fermib'
     265             : 
     266             :       INTEGER                                            :: handle, iocc, iounit, ispin, jocc, &
     267             :                                                             nactive, nspins
     268             :       REAL(KIND=dp)                                      :: dummykTS, nelec
     269          14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: interocc, occup_im
     270          14 :       REAL(kind=dp), DIMENSION(:), POINTER               :: fermia, mo_evals, occup, occuptmp
     271          14 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: fermib
     272             :       TYPE(cp_logger_type), POINTER                      :: logger
     273             :       TYPE(dft_control_type), POINTER                    :: dft_control
     274          14 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     275             :       TYPE(scf_control_type), POINTER                    :: scf_control
     276             :       TYPE(smear_type), POINTER                          :: smear
     277             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     278             : 
     279          14 :       CALL timeset(routineN, handle)
     280             : 
     281          14 :       NULLIFY (logger)
     282          14 :       logger => cp_get_default_logger()
     283          14 :       iounit = cp_logger_get_default_io_unit(logger)
     284             : 
     285          14 :       NULLIFY (mos, scf_control)
     286          14 :       CALL get_qs_env(qs_env, mos=mos, scf_control=scf_control)
     287          14 :       NULLIFY (smear)
     288          14 :       IF (ASSOCIATED(qs_env%scf_control%smear)) THEN
     289             :          smear => qs_env%scf_control%smear
     290             :       ELSE
     291           0 :          CPABORT("Smeared input section no longer associated.")
     292             :       END IF
     293             : 
     294          14 :       NULLIFY (dft_control, tddfpt_control)
     295          14 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     296          14 :       tddfpt_control => dft_control%tddfpt2_control
     297             : 
     298          14 :       NULLIFY (fermib, fermia)
     299          14 :       nspins = SIZE(gs_mos)
     300             : 
     301          28 :       DO ispin = 1, nspins
     302             : 
     303          14 :          fermib => dft_control%tddfpt2_control%smeared_occup(ispin)%fermib
     304          14 :          fermia => dft_control%tddfpt2_control%smeared_occup(ispin)%fermia
     305         294 :          fermib = 0.0_dp
     306             : 
     307             :          ! get theta_Fi
     308          14 :          NULLIFY (mo_evals, occup)
     309          14 :          CALL get_mo_set(mos(ispin), eigenvalues=mo_evals, occupation_numbers=occup)
     310          14 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nactive)
     311             : 
     312          14 :          IF (smear%fixed_mag_mom == -1.0_dp) THEN
     313           0 :             nelec = REAL(mos(ispin)%nelectron, dp)
     314             :          ELSE
     315          14 :             nelec = mos(ispin)%n_el_f
     316             :          END IF
     317             : 
     318             :          ! compute theta_im
     319          14 :          NULLIFY (occuptmp)
     320          14 :          CALL get_mo_set(mos(ispin), occupation_numbers=occuptmp)
     321          84 :          ALLOCATE (occup_im(nactive, nactive), interocc(nactive, nactive))
     322             : 
     323          70 :          DO iocc = 1, nactive
     324          56 :             IF (smear%method .EQ. smear_fermi_dirac) THEN
     325             :                ! Different prefactor in comparison to Octopus !
     326             :                CALL Fermi(occuptmp, nelec, dummykTS, mos(ispin)%eigenvalues, mos(ispin)%eigenvalues(iocc), &
     327          56 :                           smear%electronic_temperature, mos(ispin)%maxocc)
     328             :             ELSE
     329           0 :                CPABORT("TDDFT with smearing only works with Fermi-Dirac distribution.")
     330             :             END IF
     331         294 :             DO jocc = 1, nactive
     332         280 :                occup_im(iocc, jocc) = occuptmp(jocc)
     333             :             END DO
     334             :          END DO
     335             : 
     336             :          ! compute fermib
     337          70 :          DO iocc = 1, nactive
     338         294 :             DO jocc = 1, nactive
     339         224 :                interocc(iocc, jocc) = (occup(iocc) - occup(jocc))/(mo_evals(iocc) - mo_evals(jocc) - evals)
     340             :                fermib(iocc, jocc) = fermib(iocc, jocc) + occup(iocc)*occup_im(iocc, jocc) &
     341             :                                     + occup(jocc)*occup_im(jocc, iocc) &
     342         280 :                                     + fermia(jocc)*interocc(iocc, jocc)*occup_im(jocc, iocc)
     343             :             END DO
     344             :          END DO
     345             : 
     346             :          IF (debug_this_module .AND. (iounit > 0)) THEN
     347             :             DO iocc = 1, nactive
     348             :                DO jocc = 1, nactive
     349             :                   WRITE (iounit, '(A,F14.5)') "Fermi smearing parameter beta,", fermib(iocc, jocc)
     350             :                END DO
     351             :             END DO
     352             :          END IF
     353             : 
     354          14 :          DEALLOCATE (occup_im)
     355          42 :          DEALLOCATE (interocc)
     356             : 
     357             :       END DO  ! ispin
     358             : 
     359          14 :       CALL timestop(handle)
     360          28 :    END SUBROUTINE compute_fermib
     361             : ! **************************************************************************************************
     362             : !> \brief ...
     363             : !> \param qs_env ...
     364             : !> \param gs_mos ...
     365             : ! **************************************************************************************************
     366           2 :    SUBROUTINE allocate_fermi_params(qs_env, gs_mos)
     367             : 
     368             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     369             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     370             :          INTENT(IN), POINTER                             :: gs_mos
     371             : 
     372             :       CHARACTER(len=*), PARAMETER :: routineN = 'allocate_fermi_params'
     373             : 
     374             :       INTEGER                                            :: handle, ispin, nocc, nspins
     375             :       TYPE(dft_control_type), POINTER                    :: dft_control
     376           2 :       TYPE(smeared_type), DIMENSION(:), POINTER          :: smeared_occup
     377             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     378             : 
     379           2 :       CALL timeset(routineN, handle)
     380             : 
     381           2 :       NULLIFY (dft_control, tddfpt_control)
     382           2 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     383           2 :       tddfpt_control => dft_control%tddfpt2_control
     384             : 
     385             :       NULLIFY (smeared_occup)
     386           2 :       smeared_occup => dft_control%tddfpt2_control%smeared_occup
     387           2 :       nspins = SIZE(gs_mos)
     388           8 :       ALLOCATE (smeared_occup(nspins))
     389             : 
     390           4 :       DO ispin = 1, nspins
     391           2 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nocc)
     392           6 :          ALLOCATE (smeared_occup(ispin)%fermia(nocc))
     393          12 :          ALLOCATE (smeared_occup(ispin)%fermib(nocc, nocc))
     394             :       END DO
     395           2 :       dft_control%tddfpt2_control%smeared_occup => smeared_occup
     396             : 
     397           2 :       CALL timestop(handle)
     398           2 :    END SUBROUTINE allocate_fermi_params
     399             : ! **************************************************************************************************
     400             : !> \brief ...
     401             : !> \param smeared_occup ...
     402             : ! **************************************************************************************************
     403           2 :    SUBROUTINE deallocate_fermi_params(smeared_occup)
     404             : 
     405             :       TYPE(smeared_type), DIMENSION(:), POINTER          :: smeared_occup
     406             : 
     407             :       INTEGER                                            :: ispin
     408             : 
     409           2 :       IF (ASSOCIATED(smeared_occup)) THEN
     410           4 :          DO ispin = 1, SIZE(smeared_occup)
     411           4 :             IF (ASSOCIATED(smeared_occup(ispin)%fermia)) THEN
     412           2 :                DEALLOCATE (smeared_occup(ispin)%fermia)
     413           2 :                DEALLOCATE (smeared_occup(ispin)%fermib)
     414           2 :                NULLIFY (smeared_occup(ispin)%fermia, smeared_occup(ispin)%fermib)
     415             :             END IF
     416             :          END DO
     417           2 :          DEALLOCATE (smeared_occup)
     418             :          NULLIFY (smeared_occup)
     419             :       END IF
     420             : 
     421           2 :    END SUBROUTINE deallocate_fermi_params
     422             : ! **************************************************************************************************
     423             : !> \brief ...
     424             : !> \param evects ...
     425             : !> \param qs_env ...
     426             : !> \param mos_occ ...
     427             : !> \param S_C0 ...
     428             : ! **************************************************************************************************
     429          14 :    SUBROUTINE orthogonalize_smeared_occupation(evects, qs_env, mos_occ, S_C0)
     430             : 
     431             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(in)         :: evects
     432             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     433             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(in)         :: mos_occ, S_C0
     434             : 
     435             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'orthogonalize_smeared_occupation'
     436             : 
     437             :       INTEGER                                            :: handle, iocc, iounit, ispin, nactive, &
     438             :                                                             nao, nspins
     439             :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: bscale
     440          14 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: evortholocal, subevects, subevectsresult
     441          14 :       REAL(kind=dp), DIMENSION(:), POINTER               :: occup
     442          14 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: fermib
     443             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_global
     444             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     445             :       TYPE(cp_fm_type)                                   :: betaSCC, Cvec, evortho, subevectsfm, &
     446             :                                                             subevectsresultfm
     447             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     448             :       TYPE(cp_logger_type), POINTER                      :: logger
     449             :       TYPE(dft_control_type), POINTER                    :: dft_control
     450          14 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     451             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     452             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     453             : 
     454          14 :       CALL timeset(routineN, handle)
     455             : 
     456          14 :       NULLIFY (logger)
     457          14 :       logger => cp_get_default_logger()
     458          14 :       iounit = cp_logger_get_default_io_unit(logger)
     459             : 
     460          14 :       NULLIFY (mos)
     461          14 :       CALL get_qs_env(qs_env, mos=mos)
     462          14 :       NULLIFY (dft_control, tddfpt_control)
     463          14 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     464          14 :       tddfpt_control => dft_control%tddfpt2_control
     465             : 
     466             :       CALL cp_fm_get_info(matrix=evects(1), matrix_struct=matrix_struct, &
     467          14 :                           nrow_global=nao, ncol_global=nactive)
     468          14 :       CALL cp_fm_create(evortho, matrix_struct)
     469             : 
     470          14 :       nspins = SIZE(evects)
     471          14 :       NULLIFY (para_env)
     472          14 :       CALL cp_fm_get_info(evects(1), para_env=para_env, context=blacs_env_global)
     473             : 
     474          14 :       NULLIFY (matrix_struct)
     475          14 :       CALL cp_fm_struct_create(matrix_struct, nrow_global=nao, ncol_global=nao, context=blacs_env_global)
     476          14 :       CALL cp_fm_create(betaSCC, matrix_struct)
     477          14 :       CALL cp_fm_struct_release(fmstruct=matrix_struct)
     478             : 
     479          56 :       ALLOCATE (evortholocal(nao, nactive))
     480          42 :       ALLOCATE (bscale(nactive))
     481             : 
     482          28 :       DO ispin = 1, nspins
     483          14 :          NULLIFY (matrix_struct)
     484          14 :          CALL cp_fm_get_info(matrix=mos_occ(ispin), matrix_struct=matrix_struct)
     485          14 :          CALL cp_fm_create(Cvec, matrix_struct)
     486             : 
     487          14 :          NULLIFY (occup)
     488          14 :          CALL get_mo_set(mos(ispin), occupation_numbers=occup, mo_coeff=mo_coeff)
     489             : 
     490          14 :          NULLIFY (fermib)
     491          14 :          IF (.NOT. ASSOCIATED(dft_control%tddfpt2_control%smeared_occup)) THEN
     492           0 :             CPABORT("Smeared occupation intermediates not associated.")
     493             :          END IF
     494          14 :          fermib => dft_control%tddfpt2_control%smeared_occup(ispin)%fermib
     495             : 
     496          98 :          DO iocc = 1, nactive
     497          56 :             CALL cp_fm_copy_general(mos_occ(ispin), Cvec, para_env)
     498         280 :             bscale(:) = fermib(iocc, :)
     499          56 :             CALL cp_fm_column_scale(Cvec, bscale)
     500             : 
     501          56 :             CALL parallel_gemm('N', 'T', nao, nao, nactive, 1.0_dp, Cvec, S_C0(ispin), 0.0_dp, betaSCC)
     502             : 
     503             :             ! get ith column of X
     504          56 :             NULLIFY (matrix_struct)
     505          56 :             CALL cp_fm_struct_create(matrix_struct, nrow_global=nao, ncol_global=1, context=blacs_env_global)
     506          56 :             CALL cp_fm_create(subevectsfm, matrix_struct)
     507          56 :             CALL cp_fm_create(subevectsresultfm, matrix_struct)
     508          56 :             CALL cp_fm_struct_release(fmstruct=matrix_struct)
     509             : 
     510         168 :             ALLOCATE (subevects(nao, 1))
     511         112 :             ALLOCATE (subevectsresult(nao, 1))
     512             :             CALL cp_fm_get_submatrix(fm=evects(1), target_m=subevects, &
     513          56 :                                      start_row=1, start_col=iocc, n_rows=nao, n_cols=1)
     514             :             CALL cp_fm_set_submatrix(subevectsfm, subevects, &
     515          56 :                                      start_row=1, start_col=1, n_rows=nao, n_cols=1)
     516             : 
     517             :             CALL parallel_gemm('N', 'N', nao, 1, nao, 1.0_dp, betaSCC, &
     518          56 :                                subevectsfm, 0.0_dp, subevectsresultfm)
     519             : 
     520             :             CALL cp_fm_get_submatrix(fm=subevectsresultfm, target_m=subevectsresult, &
     521          56 :                                      start_row=1, start_col=1, n_rows=nao, n_cols=1)
     522             :             CALL cp_fm_set_submatrix(evortho, subevectsresult, &
     523          56 :                                      start_row=1, start_col=iocc, n_rows=nao, n_cols=1)
     524          56 :             DEALLOCATE (subevects, subevectsresult)
     525          56 :             CALL cp_fm_release(subevectsfm)
     526         126 :             CALL cp_fm_release(subevectsresultfm)
     527             :          END DO ! iocc
     528             :       END DO ! nspins
     529             : 
     530          14 :       CALL cp_fm_column_scale(evects(1), occup)
     531          14 :       CALL cp_fm_scale_and_add(1.0_dp, evects(1), -1.0_dp, evortho)
     532             : 
     533          14 :       DEALLOCATE (bscale)
     534          14 :       DEALLOCATE (evortholocal)
     535             : 
     536          14 :       CALL cp_fm_release(betaSCC)
     537          14 :       CALL cp_fm_release(Cvec)
     538             : 
     539          14 :       CALL cp_fm_release(evortho)
     540             : 
     541          14 :       CALL timestop(handle)
     542          56 :    END SUBROUTINE orthogonalize_smeared_occupation
     543             : 
     544             : END MODULE qs_tddfpt2_smearing_methods