LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_smearing_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.1 % 205 199
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            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
        

Generated by: LCOV version 2.0-1