LCOV - code coverage report
Current view: top level - src - qs_scf_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:5064cfc) Lines: 71.7 % 265 190
Test Date: 2026-03-04 06:45:10 Functions: 81.8 % 11 9

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief groups fairly general SCF methods, so that modules other than qs_scf can use them too
      10              : !>        split off from qs_scf to reduce dependencies
      11              : !> \par History
      12              : !>      - Joost VandeVondele (03.2006)
      13              : !>      - combine_ks_matrices added (05.04.06,MK)
      14              : !>      - second ROKS scheme added (15.04.06,MK)
      15              : !>      - MO occupation management moved (29.08.2008,MK)
      16              : !>      - correct_mo_eigenvalues was moved from qs_mo_types;
      17              : !>        new subroutine shift_unocc_mos (03.2016, Sergey Chulkov)
      18              : ! **************************************************************************************************
      19              : MODULE qs_scf_methods
      20              :    USE cp_dbcsr_api,                    ONLY: &
      21              :         dbcsr_add, dbcsr_desymmetrize, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
      22              :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      23              :         dbcsr_multiply, dbcsr_p_type, dbcsr_type
      24              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      25              :                                               cp_dbcsr_sm_fm_multiply
      26              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      27              :                                               cp_fm_symm,&
      28              :                                               cp_fm_triangular_multiply,&
      29              :                                               cp_fm_uplo_to_full
      30              :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_reduce,&
      31              :                                               cp_fm_cholesky_restore
      32              :    USE cp_fm_cusolver_api,              ONLY: cp_fm_general_cusolver
      33              :    USE cp_fm_diag,                      ONLY: choose_eigv_solver,&
      34              :                                               cp_fm_block_jacobi
      35              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      36              :                                               cp_fm_struct_equivalent,&
      37              :                                               cp_fm_struct_release,&
      38              :                                               cp_fm_struct_type
      39              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      40              :                                               cp_fm_get_info,&
      41              :                                               cp_fm_release,&
      42              :                                               cp_fm_to_fm,&
      43              :                                               cp_fm_type
      44              :    USE input_constants,                 ONLY: cholesky_inverse,&
      45              :                                               cholesky_off,&
      46              :                                               cholesky_reduce,&
      47              :                                               cholesky_restore
      48              :    USE kinds,                           ONLY: dp
      49              :    USE message_passing,                 ONLY: mp_para_env_type
      50              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      51              :    USE qs_density_mixing_types,         ONLY: mixing_storage_type
      52              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      53              :                                               mo_set_type
      54              : #include "./base/base_uses.f90"
      55              : 
      56              :    IMPLICIT NONE
      57              : 
      58              :    PRIVATE
      59              : 
      60              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_methods'
      61              :    REAL(KIND=dp), PARAMETER    :: ratio = 0.25_dp
      62              : 
      63              :    PUBLIC :: combine_ks_matrices, &
      64              :              cp_sm_mix, &
      65              :              eigensolver, &
      66              :              eigensolver_generalized, &
      67              :              eigensolver_dbcsr, &
      68              :              eigensolver_symm, &
      69              :              eigensolver_simple, &
      70              :              scf_env_density_mixing
      71              : 
      72              :    INTERFACE combine_ks_matrices
      73              :       MODULE PROCEDURE combine_ks_matrices_1, &
      74              :          combine_ks_matrices_2
      75              :    END INTERFACE combine_ks_matrices
      76              : 
      77              : CONTAINS
      78              : 
      79              : ! **************************************************************************************************
      80              : !> \brief perform (if requested) a density mixing
      81              : !> \param p_mix_new    New density matrices
      82              : !> \param mixing_store ...
      83              : !> \param rho_ao       Density environment
      84              : !> \param para_env ...
      85              : !> \param iter_delta ...
      86              : !> \param iter_count ...
      87              : !> \param diis ...
      88              : !> \param invert       Invert mixing
      89              : !> \par History
      90              : !>      02.2003 created [fawzi]
      91              : !>      08.2014 adapted for kpoints [JGH]
      92              : !> \author fawzi
      93              : ! **************************************************************************************************
      94       110395 :    SUBROUTINE scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, &
      95              :                                      iter_delta, iter_count, diis, invert)
      96              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: p_mix_new
      97              :       TYPE(mixing_storage_type), POINTER                 :: mixing_store
      98              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao
      99              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     100              :       REAL(KIND=dp), INTENT(INOUT)                       :: iter_delta
     101              :       INTEGER, INTENT(IN)                                :: iter_count
     102              :       LOGICAL, INTENT(in), OPTIONAL                      :: diis, invert
     103              : 
     104              :       CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_density_mixing'
     105              : 
     106              :       INTEGER                                            :: handle, ic, ispin
     107              :       LOGICAL                                            :: my_diis, my_invert
     108              :       REAL(KIND=dp)                                      :: my_p_mix, tmp
     109              : 
     110       110395 :       CALL timeset(routineN, handle)
     111              : 
     112       110395 :       my_diis = .FALSE.
     113       110395 :       IF (PRESENT(diis)) my_diis = diis
     114       110395 :       my_invert = .FALSE.
     115       110395 :       IF (PRESENT(invert)) my_invert = invert
     116       110395 :       my_p_mix = mixing_store%alpha
     117       110395 :       IF (my_diis .OR. iter_count < mixing_store%nskip_mixing) THEN
     118        71781 :          my_p_mix = 1.0_dp
     119              :       END IF
     120              : 
     121       110395 :       iter_delta = 0.0_dp
     122       110395 :       CPASSERT(ASSOCIATED(p_mix_new))
     123       534392 :       DO ic = 1, SIZE(p_mix_new, 2)
     124      1032379 :          DO ispin = 1, SIZE(p_mix_new, 1)
     125       921984 :             IF (my_invert) THEN
     126        68352 :                CPASSERT(my_p_mix /= 0.0_dp)
     127        68352 :                IF (my_p_mix /= 1.0_dp) THEN
     128              :                   CALL dbcsr_add(matrix_a=p_mix_new(ispin, ic)%matrix, &
     129              :                                  alpha_scalar=1.0_dp/my_p_mix, &
     130              :                                  matrix_b=rho_ao(ispin, ic)%matrix, &
     131        16008 :                                  beta_scalar=(my_p_mix - 1.0_dp)/my_p_mix)
     132              :                END IF
     133              :             ELSE
     134              :                CALL cp_sm_mix(m1=p_mix_new(ispin, ic)%matrix, &
     135              :                               m2=rho_ao(ispin, ic)%matrix, &
     136              :                               p_mix=my_p_mix, &
     137              :                               delta=tmp, &
     138       429635 :                               para_env=para_env)
     139       429635 :                iter_delta = MAX(iter_delta, tmp)
     140              :             END IF
     141              :          END DO
     142              :       END DO
     143              : 
     144       110395 :       CALL timestop(handle)
     145              : 
     146       110395 :    END SUBROUTINE scf_env_density_mixing
     147              : 
     148              : ! **************************************************************************************************
     149              : !> \brief   Diagonalise the Kohn-Sham matrix to get a new set of MO eigen-
     150              : !>          vectors and MO eigenvalues. ks will be modified
     151              : !> \param matrix_ks_fm ...
     152              : !> \param mo_set ...
     153              : !> \param ortho ...
     154              : !> \param work ...
     155              : !> \param cholesky_method ...
     156              : !> \param do_level_shift activate the level shifting technique
     157              : !> \param level_shift    amount of shift applied (in a.u.)
     158              : !> \param matrix_u_fm    matrix U : S (overlap matrix) = U^T * U
     159              : !> \param use_jacobi ...
     160              : !> \date    01.05.2001
     161              : !> \author  Matthias Krack
     162              : !> \version 1.0
     163              : ! **************************************************************************************************
     164       173918 :    SUBROUTINE eigensolver(matrix_ks_fm, mo_set, ortho, work, &
     165              :                           cholesky_method, do_level_shift, &
     166              :                           level_shift, matrix_u_fm, use_jacobi)
     167              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_ks_fm
     168              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     169              :       TYPE(cp_fm_type), INTENT(IN)                       :: ortho, work
     170              :       INTEGER, INTENT(inout)                             :: cholesky_method
     171              :       LOGICAL, INTENT(in)                                :: do_level_shift
     172              :       REAL(KIND=dp), INTENT(in)                          :: level_shift
     173              :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: matrix_u_fm
     174              :       LOGICAL, INTENT(in)                                :: use_jacobi
     175              : 
     176              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'eigensolver'
     177              : 
     178              :       INTEGER                                            :: handle, homo, nao, nmo
     179        86959 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues
     180              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     181              : 
     182        86959 :       CALL timeset(routineN, handle)
     183              : 
     184        86959 :       NULLIFY (mo_coeff)
     185        86959 :       NULLIFY (mo_eigenvalues)
     186              : 
     187              :       ! Diagonalise the Kohn-Sham matrix
     188              : 
     189              :       CALL get_mo_set(mo_set=mo_set, &
     190              :                       nao=nao, &
     191              :                       nmo=nmo, &
     192              :                       homo=homo, &
     193              :                       eigenvalues=mo_eigenvalues, &
     194        86959 :                       mo_coeff=mo_coeff)
     195              : 
     196        86999 :       SELECT CASE (cholesky_method)
     197              :       CASE (cholesky_reduce)
     198           40 :          CALL cp_fm_cholesky_reduce(matrix_ks_fm, ortho)
     199              : 
     200           40 :          IF (do_level_shift) &
     201              :             CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
     202           28 :                                  level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
     203              : 
     204           40 :          CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
     205           40 :          CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
     206           40 :          IF (do_level_shift) &
     207           28 :             CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
     208              : 
     209              :       CASE (cholesky_restore)
     210        85941 :          CALL cp_fm_uplo_to_full(matrix_ks_fm, work)
     211              :          CALL cp_fm_cholesky_restore(matrix_ks_fm, nao, ortho, work, &
     212        85941 :                                      "SOLVE", pos="RIGHT")
     213              :          CALL cp_fm_cholesky_restore(work, nao, ortho, matrix_ks_fm, &
     214        85941 :                                      "SOLVE", pos="LEFT", transa="T")
     215              : 
     216        85941 :          IF (do_level_shift) &
     217              :             CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
     218           68 :                                  level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
     219              : 
     220        85941 :          CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
     221        85941 :          CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
     222              : 
     223        85941 :          IF (do_level_shift) &
     224           68 :             CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
     225              : 
     226              :       CASE (cholesky_inverse)
     227          978 :          CALL cp_fm_uplo_to_full(matrix_ks_fm, work)
     228              : 
     229              :          CALL cp_fm_triangular_multiply(ortho, matrix_ks_fm, side="R", transpose_tr=.FALSE., &
     230          978 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
     231              :          CALL cp_fm_triangular_multiply(ortho, matrix_ks_fm, side="L", transpose_tr=.TRUE., &
     232          978 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
     233              : 
     234          978 :          IF (do_level_shift) &
     235              :             CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
     236           28 :                                  level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
     237              : 
     238          978 :          CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
     239              :          CALL cp_fm_triangular_multiply(ortho, work, side="L", transpose_tr=.FALSE., &
     240          978 :                                         invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
     241          978 :          CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
     242              : 
     243          978 :          IF (do_level_shift) &
     244        86987 :             CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
     245              : 
     246              :       END SELECT
     247              : 
     248        86959 :       IF (use_jacobi) THEN
     249            0 :          CALL cp_fm_to_fm(mo_coeff, ortho)
     250            0 :          cholesky_method = cholesky_off
     251              :       END IF
     252              : 
     253        86959 :       CALL timestop(handle)
     254              : 
     255        86959 :    END SUBROUTINE eigensolver
     256              : 
     257              : ! **************************************************************************************************
     258              : !> \brief Solve the generalized eigenvalue problem using cusolverMpSygvd
     259              : !> \param matrix_ks_fm Kohn-Sham matrix in FM format
     260              : !> \param matrix_s     Overlap matrix (DBCSR)
     261              : !> \param mo_set       Molecular orbital set
     262              : !> \param work         Work matrix (used as eigenvector buffer)
     263              : ! **************************************************************************************************
     264            0 :    SUBROUTINE eigensolver_generalized(matrix_ks_fm, matrix_s, mo_set, work)
     265              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: matrix_ks_fm
     266              :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_s
     267              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     268              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: work
     269              : 
     270              :       CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_generalized'
     271              : 
     272              :       INTEGER                                            :: handle, nmo
     273            0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues
     274              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     275              :       TYPE(cp_fm_type)                                   :: s_fm
     276              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     277              : 
     278            0 :       CALL timeset(routineN, handle)
     279              : 
     280            0 :       NULLIFY (mo_coeff)
     281            0 :       NULLIFY (mo_eigenvalues)
     282              : 
     283            0 :       CALL get_mo_set(mo_set=mo_set, nmo=nmo, eigenvalues=mo_eigenvalues, mo_coeff=mo_coeff)
     284            0 :       CALL cp_fm_get_info(matrix_ks_fm, matrix_struct=fm_struct)
     285              : 
     286              :       ! Convert S matrix from DBCSR to FM (required for cuSOLVERMp)
     287            0 :       CALL cp_fm_create(s_fm, fm_struct)
     288            0 :       CALL copy_dbcsr_to_fm(matrix_s, s_fm)
     289              : 
     290              :       ! Solve generalized eigenvalue problem - eigenvectors output to work buffer
     291            0 :       CALL cp_fm_general_cusolver(matrix_ks_fm, s_fm, work, mo_eigenvalues)
     292              : 
     293              :       ! Copy only the occupied MOs to mo_coeff
     294            0 :       CALL cp_fm_to_fm(work, mo_coeff, nmo)
     295              : 
     296            0 :       CALL cp_fm_release(s_fm)
     297              : 
     298            0 :       CALL timestop(handle)
     299              : 
     300            0 :    END SUBROUTINE eigensolver_generalized
     301              : 
     302              : ! **************************************************************************************************
     303              : !> \brief ...
     304              : !> \param matrix_ks ...
     305              : !> \param matrix_ks_fm ...
     306              : !> \param mo_set ...
     307              : !> \param ortho_dbcsr ...
     308              : !> \param ksbuf1 ...
     309              : !> \param ksbuf2 ...
     310              : ! **************************************************************************************************
     311         8472 :    SUBROUTINE eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
     312              :       TYPE(dbcsr_type), INTENT(IN)                       :: matrix_ks
     313              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: matrix_ks_fm
     314              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     315              :       TYPE(dbcsr_type), INTENT(IN)                       :: ortho_dbcsr
     316              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: ksbuf1, ksbuf2
     317              : 
     318              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'eigensolver_dbcsr'
     319              : 
     320              :       INTEGER                                            :: handle, nao, nmo
     321         2118 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues
     322              :       TYPE(cp_fm_type)                                   :: all_evecs, nmo_evecs
     323              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     324              : 
     325         2118 :       CALL timeset(routineN, handle)
     326              : 
     327         2118 :       NULLIFY (mo_coeff)
     328         2118 :       NULLIFY (mo_eigenvalues)
     329              : 
     330              :       CALL get_mo_set(mo_set=mo_set, &
     331              :                       nao=nao, &
     332              :                       nmo=nmo, &
     333              :                       eigenvalues=mo_eigenvalues, &
     334         2118 :                       mo_coeff=mo_coeff)
     335              : 
     336              : !    Reduce KS matrix
     337         2118 :       CALL dbcsr_desymmetrize(matrix_ks, ksbuf2)
     338         2118 :       CALL dbcsr_multiply('N', 'N', 1.0_dp, ksbuf2, ortho_dbcsr, 0.0_dp, ksbuf1)
     339         2118 :       CALL dbcsr_multiply('T', 'N', 1.0_dp, ortho_dbcsr, ksbuf1, 0.0_dp, ksbuf2)
     340              : 
     341              : !    Solve the eigenvalue problem
     342         2118 :       CALL copy_dbcsr_to_fm(ksbuf2, matrix_ks_fm)
     343         2118 :       CALL cp_fm_create(all_evecs, matrix_ks_fm%matrix_struct)
     344         2118 :       CALL choose_eigv_solver(matrix_ks_fm, all_evecs, mo_eigenvalues)
     345              : 
     346              :       ! select first nmo eigenvectors
     347         2118 :       CALL cp_fm_create(nmo_evecs, mo_coeff%matrix_struct)
     348         2118 :       CALL cp_fm_to_fm(msource=all_evecs, mtarget=nmo_evecs, ncol=nmo)
     349         2118 :       CALL cp_fm_release(all_evecs)
     350              : 
     351              : !    Restore the eigenvector of the general eig. problem
     352         2118 :       CALL cp_dbcsr_sm_fm_multiply(ortho_dbcsr, nmo_evecs, mo_coeff, nmo)
     353              : 
     354         2118 :       CALL cp_fm_release(nmo_evecs)
     355         2118 :       CALL timestop(handle)
     356              : 
     357         2118 :    END SUBROUTINE eigensolver_dbcsr
     358              : 
     359              : ! **************************************************************************************************
     360              : !> \brief ...
     361              : !> \param matrix_ks_fm ...
     362              : !> \param mo_set ...
     363              : !> \param ortho ...
     364              : !> \param work ...
     365              : !> \param do_level_shift activate the level shifting technique
     366              : !> \param level_shift    amount of shift applied (in a.u.)
     367              : !> \param matrix_u_fm    matrix U : S (overlap matrix) = U^T * U
     368              : !> \param use_jacobi ...
     369              : !> \param jacobi_threshold ...
     370              : !> \param ortho_red ...
     371              : !> \param work_red ...
     372              : !> \param matrix_ks_fm_red ...
     373              : !> \param matrix_u_fm_red ...
     374              : ! **************************************************************************************************
     375          524 :    SUBROUTINE eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, &
     376              :                                level_shift, matrix_u_fm, use_jacobi, jacobi_threshold, &
     377              :                                ortho_red, work_red, matrix_ks_fm_red, matrix_u_fm_red)
     378              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_ks_fm
     379              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     380              :       TYPE(cp_fm_type), INTENT(IN)                       :: ortho, work
     381              :       LOGICAL, INTENT(IN)                                :: do_level_shift
     382              :       REAL(KIND=dp), INTENT(IN)                          :: level_shift
     383              :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: matrix_u_fm
     384              :       LOGICAL, INTENT(IN)                                :: use_jacobi
     385              :       REAL(KIND=dp), INTENT(IN)                          :: jacobi_threshold
     386              :       TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL          :: ortho_red, work_red, matrix_ks_fm_red, &
     387              :                                                             matrix_u_fm_red
     388              : 
     389              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'eigensolver_symm'
     390              : 
     391              :       INTEGER                                            :: handle, homo, nao, nao_red, nelectron, &
     392              :                                                             nmo
     393          524 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues
     394          524 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues
     395              :       TYPE(cp_fm_type)                                   :: work_red2
     396              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     397              : 
     398          524 :       CALL timeset(routineN, handle)
     399              : 
     400          524 :       NULLIFY (mo_coeff)
     401          524 :       NULLIFY (mo_eigenvalues)
     402              : 
     403              :       ! Diagonalise the Kohn-Sham matrix
     404              : 
     405              :       CALL get_mo_set(mo_set=mo_set, &
     406              :                       nao=nao, &
     407              :                       nmo=nmo, &
     408              :                       homo=homo, &
     409              :                       nelectron=nelectron, &
     410              :                       eigenvalues=mo_eigenvalues, &
     411          524 :                       mo_coeff=mo_coeff)
     412              : 
     413          524 :       IF (use_jacobi) THEN
     414              : 
     415            0 :          CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks_fm, mo_coeff, 0.0_dp, work)
     416              :          CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
     417            0 :                             0.0_dp, matrix_ks_fm, b_first_col=homo + 1, c_first_col=homo + 1)
     418              : 
     419              :          ! Block Jacobi (pseudo-diagonalization, only one sweep)
     420              :          CALL cp_fm_block_jacobi(matrix_ks_fm, mo_coeff, mo_eigenvalues, &
     421            0 :                                  jacobi_threshold, homo + 1)
     422              : 
     423              :       ELSE ! full S^(-1/2) has been computed
     424          524 :          IF (PRESENT(work_red) .AND. PRESENT(ortho_red) .AND. PRESENT(matrix_ks_fm_red)) THEN
     425          524 :             CALL cp_fm_get_info(ortho_red, ncol_global=nao_red)
     426          524 :             CALL cp_fm_symm("L", "U", nao, nao_red, 1.0_dp, matrix_ks_fm, ortho_red, 0.0_dp, work_red)
     427          524 :             CALL parallel_gemm("T", "N", nao_red, nao_red, nao, 1.0_dp, ortho_red, work_red, 0.0_dp, matrix_ks_fm_red)
     428              : 
     429          524 :             IF (do_level_shift) &
     430              :                CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm_red, mo_coeff=mo_coeff, homo=homo, &
     431           86 :                                     level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm_red)
     432              : 
     433          524 :             CALL cp_fm_create(work_red2, matrix_ks_fm_red%matrix_struct)
     434         1572 :             ALLOCATE (eigenvalues(nao_red))
     435          524 :             CALL choose_eigv_solver(matrix_ks_fm_red, work_red2, eigenvalues)
     436         6812 :             mo_eigenvalues(1:MIN(nao_red, nmo)) = eigenvalues(1:MIN(nao_red, nmo))
     437              :             CALL parallel_gemm("N", "N", nao, nmo, nao_red, 1.0_dp, ortho_red, work_red2, 0.0_dp, &
     438          524 :                                mo_coeff)
     439         1572 :             CALL cp_fm_release(work_red2)
     440              :          ELSE
     441            0 :             CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, matrix_ks_fm, ortho, 0.0_dp, work)
     442            0 :             CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, matrix_ks_fm)
     443            0 :             IF (do_level_shift) &
     444              :                CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
     445            0 :                                     level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm)
     446            0 :             CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
     447              :             CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, &
     448            0 :                                mo_coeff)
     449              :          END IF
     450              : 
     451          524 :          IF (do_level_shift) &
     452           86 :             CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
     453              : 
     454              :       END IF
     455              : 
     456          524 :       CALL timestop(handle)
     457              : 
     458         1048 :    END SUBROUTINE eigensolver_symm
     459              : 
     460              : ! **************************************************************************************************
     461              : 
     462              : ! **************************************************************************************************
     463              : !> \brief ...
     464              : !> \param matrix_ks ...
     465              : !> \param mo_set ...
     466              : !> \param work ...
     467              : !> \param do_level_shift activate the level shifting technique
     468              : !> \param level_shift    amount of shift applied (in a.u.)
     469              : !> \param use_jacobi ...
     470              : !> \param jacobi_threshold ...
     471              : ! **************************************************************************************************
     472        37356 :    SUBROUTINE eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, &
     473              :                                  level_shift, use_jacobi, jacobi_threshold)
     474              : 
     475              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_ks
     476              :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     477              :       TYPE(cp_fm_type), INTENT(IN)                       :: work
     478              :       LOGICAL, INTENT(IN)                                :: do_level_shift
     479              :       REAL(KIND=dp), INTENT(IN)                          :: level_shift
     480              :       LOGICAL, INTENT(IN)                                :: use_jacobi
     481              :       REAL(KIND=dp), INTENT(IN)                          :: jacobi_threshold
     482              : 
     483              :       CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_simple'
     484              : 
     485              :       INTEGER                                            :: handle, homo, nao, nelectron, nmo
     486        18678 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mo_eigenvalues
     487              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     488              : 
     489        18678 :       CALL timeset(routineN, handle)
     490              : 
     491        18678 :       NULLIFY (mo_coeff)
     492        18678 :       NULLIFY (mo_eigenvalues)
     493              : 
     494              :       CALL get_mo_set(mo_set=mo_set, &
     495              :                       nao=nao, &
     496              :                       nmo=nmo, &
     497              :                       homo=homo, &
     498              :                       nelectron=nelectron, &
     499              :                       eigenvalues=mo_eigenvalues, &
     500        18678 :                       mo_coeff=mo_coeff)
     501              : 
     502        18678 :       IF (do_level_shift) THEN
     503              :          ! matrix_u_fm is simply an identity matrix, so we omit it here
     504              :          CALL shift_unocc_mos(matrix_ks_fm=matrix_ks, mo_coeff=mo_coeff, homo=homo, &
     505            0 :                               level_shift=level_shift, is_triangular=.FALSE.)
     506              :       END IF
     507              : 
     508        18678 :       IF (use_jacobi) THEN
     509           18 :          CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks, mo_coeff, 0.0_dp, work)
     510              :          CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
     511           18 :                             0.0_dp, matrix_ks, b_first_col=homo + 1, c_first_col=homo + 1)
     512              :          ! Block Jacobi (pseudo-diagonalization, only one sweep)
     513           18 :          CALL cp_fm_block_jacobi(matrix_ks, mo_coeff, mo_eigenvalues, jacobi_threshold, homo + 1)
     514              :       ELSE
     515              : 
     516        18660 :          CALL choose_eigv_solver(matrix_ks, work, mo_eigenvalues)
     517              : 
     518        18660 :          CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
     519              : 
     520              :       END IF
     521              : 
     522        18678 :       IF (do_level_shift) &
     523            0 :          CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
     524              : 
     525        18678 :       CALL timestop(handle)
     526              : 
     527        18678 :    END SUBROUTINE eigensolver_simple
     528              : 
     529              : ! **************************************************************************************************
     530              : !> \brief Perform a mixing of the given matrixes into the first matrix
     531              : !>      m1 = m2 + p_mix (m1-m2)
     532              : !> \param m1 first (new) matrix, is modified
     533              : !> \param m2 the second (old) matrix
     534              : !> \param p_mix how much m1 is conserved (0: none, 1: all)
     535              : !> \param delta maximum norm of m1-m2
     536              : !> \param para_env ...
     537              : !> \param m3 ...
     538              : !> \par History
     539              : !>      02.2003 rewamped [fawzi]
     540              : !> \author fawzi
     541              : !> \note
     542              : !>      if you what to store the result in m2 swap m1 and m2 an use
     543              : !>      (1-pmix) as pmix
     544              : !>      para_env should be removed (embedded in matrix)
     545              : ! **************************************************************************************************
     546      1160030 :    SUBROUTINE cp_sm_mix(m1, m2, p_mix, delta, para_env, m3)
     547              : 
     548              :       TYPE(dbcsr_type), POINTER                          :: m1, m2
     549              :       REAL(KIND=dp), INTENT(IN)                          :: p_mix
     550              :       REAL(KIND=dp), INTENT(OUT)                         :: delta
     551              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     552              :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: m3
     553              : 
     554              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_sm_mix'
     555              : 
     556              :       INTEGER                                            :: handle, i, iblock_col, iblock_row, j
     557              :       LOGICAL                                            :: found
     558       580015 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_delta_block, p_new_block, p_old_block
     559              :       TYPE(dbcsr_iterator_type)                          :: iter
     560              : 
     561       580015 :       CALL timeset(routineN, handle)
     562       580015 :       delta = 0.0_dp
     563              : 
     564       580015 :       CALL dbcsr_iterator_start(iter, m1)
     565     12594582 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     566     12014567 :          CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_new_block)
     567              :          CALL dbcsr_get_block_p(matrix=m2, row=iblock_row, col=iblock_col, &
     568     12014567 :                                 BLOCK=p_old_block, found=found)
     569     12014567 :          CPASSERT(ASSOCIATED(p_old_block))
     570     12594582 :          IF (PRESENT(m3)) THEN
     571              :             CALL dbcsr_get_block_p(matrix=m3, row=iblock_row, col=iblock_col, &
     572      2474503 :                                    BLOCK=p_delta_block, found=found)
     573      2474503 :             CPASSERT(ASSOCIATED(p_delta_block))
     574              : 
     575     24127928 :             DO j = 1, SIZE(p_new_block, 2)
     576    221977223 :                DO i = 1, SIZE(p_new_block, 1)
     577    197849295 :                   p_delta_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
     578    219502720 :                   delta = MAX(delta, ABS(p_delta_block(i, j)))
     579              :                END DO
     580              :             END DO
     581              :          ELSE
     582     42561437 :             DO j = 1, SIZE(p_new_block, 2)
     583    267968728 :                DO i = 1, SIZE(p_new_block, 1)
     584    225407291 :                   p_new_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
     585    225407291 :                   delta = MAX(delta, ABS(p_new_block(i, j)))
     586    258428664 :                   p_new_block(i, j) = p_old_block(i, j) + p_mix*p_new_block(i, j)
     587              :                END DO
     588              :             END DO
     589              :          END IF
     590              :       END DO
     591       580015 :       CALL dbcsr_iterator_stop(iter)
     592              : 
     593       580015 :       CALL para_env%max(delta)
     594              : 
     595       580015 :       CALL timestop(handle)
     596              : 
     597       580015 :    END SUBROUTINE cp_sm_mix
     598              : 
     599              : ! **************************************************************************************************
     600              : !> \brief ...
     601              : !> \param ksa ...
     602              : !> \param ksb ...
     603              : !> \param occa ...
     604              : !> \param occb ...
     605              : !> \param roks_parameter ...
     606              : ! **************************************************************************************************
     607          998 :    SUBROUTINE combine_ks_matrices_1(ksa, ksb, occa, occb, roks_parameter)
     608              : 
     609              :       ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
     610              :       ! Kohn-Sham (ROKS) calculation
     611              :       ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
     612              :       ! respectively. occa and occb contain the corresponding MO occupation
     613              :       ! numbers. On output the combined ROKS operator matrix is returned in ksa.
     614              : 
     615              :       ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
     616              :       !             - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
     617              : 
     618              :       TYPE(cp_fm_type), INTENT(IN)                       :: ksa, ksb
     619              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: occa, occb
     620              :       REAL(KIND=dp), DIMENSION(0:2, 0:2, 1:2), &
     621              :          INTENT(IN)                                      :: roks_parameter
     622              : 
     623              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_1'
     624              : 
     625              :       INTEGER                                            :: handle, i, icol_global, icol_local, &
     626              :                                                             irow_global, irow_local, j, &
     627              :                                                             ncol_local, nrow_local
     628          998 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     629              :       LOGICAL                                            :: compatible_matrices
     630              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
     631          998 :          POINTER                                         :: fa, fb
     632              :       TYPE(cp_fm_struct_type), POINTER                   :: ksa_struct, ksb_struct
     633              : 
     634              : ! -------------------------------------------------------------------------
     635              : 
     636          998 :       CALL timeset(routineN, handle)
     637              : 
     638              :       CALL cp_fm_get_info(matrix=ksa, &
     639              :                           matrix_struct=ksa_struct, &
     640              :                           nrow_local=nrow_local, &
     641              :                           ncol_local=ncol_local, &
     642              :                           row_indices=row_indices, &
     643              :                           col_indices=col_indices, &
     644          998 :                           local_data=fa)
     645              : 
     646              :       CALL cp_fm_get_info(matrix=ksb, &
     647              :                           matrix_struct=ksb_struct, &
     648          998 :                           local_data=fb)
     649              : 
     650          998 :       compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
     651          998 :       CPASSERT(compatible_matrices)
     652              : 
     653        42085 :       IF (SUM(occb) == 0.0_dp) fb = 0.0_dp
     654              : 
     655        16196 :       DO icol_local = 1, ncol_local
     656        15198 :          icol_global = col_indices(icol_local)
     657        15198 :          j = INT(occa(icol_global)) + INT(occb(icol_global))
     658       180579 :          DO irow_local = 1, nrow_local
     659       164383 :             irow_global = row_indices(irow_local)
     660       164383 :             i = INT(occa(irow_global)) + INT(occb(irow_global))
     661              :             fa(irow_local, icol_local) = &
     662              :                roks_parameter(i, j, 1)*fa(irow_local, icol_local) + &
     663       179581 :                roks_parameter(i, j, 2)*fb(irow_local, icol_local)
     664              :          END DO
     665              :       END DO
     666              : 
     667          998 :       CALL timestop(handle)
     668              : 
     669          998 :    END SUBROUTINE combine_ks_matrices_1
     670              : 
     671              : ! **************************************************************************************************
     672              : !> \brief ...
     673              : !> \param ksa ...
     674              : !> \param ksb ...
     675              : !> \param occa ...
     676              : !> \param occb ...
     677              : !> \param f ...
     678              : !> \param nalpha ...
     679              : !> \param nbeta ...
     680              : ! **************************************************************************************************
     681            0 :    SUBROUTINE combine_ks_matrices_2(ksa, ksb, occa, occb, f, nalpha, nbeta)
     682              : 
     683              :       ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
     684              :       ! Kohn-Sham (ROKS) calculation
     685              :       ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
     686              :       ! respectively. occa and occb contain the corresponding MO occupation
     687              :       ! numbers. On output the combined ROKS operator matrix is returned in ksa.
     688              : 
     689              :       ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
     690              :       !             - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
     691              : 
     692              :       TYPE(cp_fm_type), INTENT(IN)                       :: ksa, ksb
     693              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: occa, occb
     694              :       REAL(KIND=dp), INTENT(IN)                          :: f
     695              :       INTEGER, INTENT(IN)                                :: nalpha, nbeta
     696              : 
     697              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_2'
     698              : 
     699              :       INTEGER                                            :: handle, icol_global, icol_local, &
     700              :                                                             irow_global, irow_local, ncol_local, &
     701              :                                                             nrow_local
     702            0 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     703              :       LOGICAL                                            :: compatible_matrices
     704              :       REAL(KIND=dp)                                      :: beta, t1, t2, ta, tb
     705              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
     706            0 :          POINTER                                         :: fa, fb
     707              :       TYPE(cp_fm_struct_type), POINTER                   :: ksa_struct, ksb_struct
     708              : 
     709              : ! -------------------------------------------------------------------------
     710              : 
     711            0 :       CALL timeset(routineN, handle)
     712              : 
     713              :       CALL cp_fm_get_info(matrix=ksa, &
     714              :                           matrix_struct=ksa_struct, &
     715              :                           nrow_local=nrow_local, &
     716              :                           ncol_local=ncol_local, &
     717              :                           row_indices=row_indices, &
     718              :                           col_indices=col_indices, &
     719            0 :                           local_data=fa)
     720              : 
     721              :       CALL cp_fm_get_info(matrix=ksb, &
     722              :                           matrix_struct=ksb_struct, &
     723            0 :                           local_data=fb)
     724              : 
     725            0 :       compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
     726            0 :       CPASSERT(compatible_matrices)
     727              : 
     728            0 :       beta = 1.0_dp/(1.0_dp - f)
     729              : 
     730            0 :       DO icol_local = 1, ncol_local
     731              : 
     732            0 :          icol_global = col_indices(icol_local)
     733              : 
     734            0 :          DO irow_local = 1, nrow_local
     735              : 
     736            0 :             irow_global = row_indices(irow_local)
     737              : 
     738            0 :             t1 = 0.5_dp*(fa(irow_local, icol_local) + fb(irow_local, icol_local))
     739              : 
     740            0 :             IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
     741            0 :                IF ((0 < icol_global) .AND. (icol_global <= nbeta)) THEN
     742              :                   ! closed-closed
     743            0 :                   fa(irow_local, icol_local) = t1
     744            0 :                ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
     745              :                   ! closed-open
     746            0 :                   ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
     747            0 :                   tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
     748            0 :                   t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
     749            0 :                   fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
     750              :                ELSE
     751              :                   ! closed-virtual
     752            0 :                   fa(irow_local, icol_local) = t1
     753              :                END IF
     754            0 :             ELSE IF ((nbeta < irow_global) .AND. (irow_global <= nalpha)) THEN
     755              :                IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
     756              :                   ! open-closed
     757              :                   ta = 0.5_dp*(f - REAL(occa(irow_global), KIND=dp))/f
     758              :                   tb = 0.5_dp*(f - REAL(occb(irow_global), KIND=dp))/f
     759              :                   t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
     760              :                   fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
     761            0 :                ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
     762              :                   ! open-open
     763            0 :                   ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
     764            0 :                   tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
     765            0 :                   t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
     766            0 :                   IF (irow_global == icol_global) THEN
     767            0 :                      fa(irow_local, icol_local) = t1 - t2
     768              :                   ELSE
     769            0 :                      fa(irow_local, icol_local) = t1 - 0.5_dp*t2
     770              :                   END IF
     771              :                ELSE
     772              :                   ! open-virtual
     773            0 :                   ta = 0.5_dp*(f - REAL(occa(irow_global), KIND=dp))/f
     774            0 :                   tb = 0.5_dp*(f - REAL(occb(irow_global), KIND=dp))/f
     775            0 :                   t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
     776            0 :                   fa(irow_local, icol_local) = t1 - t2
     777              :                END IF
     778              :             ELSE
     779            0 :                IF ((0 < irow_global) .AND. (irow_global < nbeta)) THEN
     780              :                   ! virtual-closed
     781            0 :                   fa(irow_local, icol_local) = t1
     782            0 :                ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
     783              :                   ! virtual-open
     784            0 :                   ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
     785            0 :                   tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
     786            0 :                   t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
     787            0 :                   fa(irow_local, icol_local) = t1 - t2
     788              :                ELSE
     789              :                   ! virtual-virtual
     790            0 :                   fa(irow_local, icol_local) = t1
     791              :                END IF
     792              :             END IF
     793              : 
     794              :          END DO
     795              :       END DO
     796              : 
     797            0 :       CALL timestop(handle)
     798              : 
     799            0 :    END SUBROUTINE combine_ks_matrices_2
     800              : 
     801              : ! **************************************************************************************************
     802              : !> \brief   Correct MO eigenvalues after MO level shifting.
     803              : !> \param mo_eigenvalues vector of eigenvalues
     804              : !> \param homo index of the highest occupied molecular orbital
     805              : !> \param nmo  number of molecular orbitals
     806              : !> \param level_shift amount of applied level shifting (in a.u.)
     807              : !> \date    19.04.2002
     808              : !> \par History
     809              : !>      - correct_mo_eigenvalues added (18.04.02,MK)
     810              : !>      - moved from module qs_mo_types, revised interface (03.2016, Sergey Chulkov)
     811              : !> \author  MK
     812              : !> \version 1.0
     813              : ! **************************************************************************************************
     814          210 :    PURE SUBROUTINE correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
     815              : 
     816              :       REAL(kind=dp), DIMENSION(:), INTENT(inout)         :: mo_eigenvalues
     817              :       INTEGER, INTENT(in)                                :: homo, nmo
     818              :       REAL(kind=dp), INTENT(in)                          :: level_shift
     819              : 
     820              :       INTEGER                                            :: imo
     821              : 
     822         5530 :       DO imo = homo + 1, nmo
     823         5530 :          mo_eigenvalues(imo) = mo_eigenvalues(imo) - level_shift
     824              :       END DO
     825              : 
     826          210 :    END SUBROUTINE correct_mo_eigenvalues
     827              : 
     828              : ! **************************************************************************************************
     829              : !> \brief Adjust the Kohn-Sham matrix by shifting the orbital energies of all
     830              : !>        unoccupied molecular orbitals
     831              : !> \param matrix_ks_fm   transformed Kohn-Sham matrix = U^{-1,T} * KS * U^{-1}
     832              : !> \param mo_coeff       matrix of molecular orbitals (C)
     833              : !> \param homo           number of occupied molecular orbitals
     834              : !> \param level_shift    amount of shift applying (in a.u.)
     835              : !> \param is_triangular  indicates that matrix_u_fm contains an upper triangular matrix
     836              : !> \param matrix_u_fm    matrix U: S (overlap matrix) = U^T * U;
     837              : !>                       assume an identity matrix if omitted
     838              : !> \par History
     839              : !>      03.2016 created [Sergey Chulkov]
     840              : ! **************************************************************************************************
     841          210 :    SUBROUTINE shift_unocc_mos(matrix_ks_fm, mo_coeff, homo, &
     842              :                               level_shift, is_triangular, matrix_u_fm)
     843              : 
     844              :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_ks_fm, mo_coeff
     845              :       INTEGER, INTENT(in)                                :: homo
     846              :       REAL(kind=dp), INTENT(in)                          :: level_shift
     847              :       LOGICAL, INTENT(in)                                :: is_triangular
     848              :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: matrix_u_fm
     849              : 
     850              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'shift_unocc_mos'
     851              : 
     852              :       INTEGER                                            :: handle, nao, nao_red, nmo
     853          210 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: weights
     854              :       TYPE(cp_fm_struct_type), POINTER                   :: ao_mo_fmstruct
     855              :       TYPE(cp_fm_type)                                   :: u_mo, u_mo_scaled
     856              : 
     857          210 :       CALL timeset(routineN, handle)
     858              : 
     859          210 :       IF (PRESENT(matrix_u_fm)) THEN
     860          210 :          CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
     861          210 :          CALL cp_fm_get_info(matrix_u_fm, nrow_global=nao_red, ncol_global=nao)
     862              :       ELSE
     863            0 :          CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
     864            0 :          nao_red = nao
     865              :       END IF
     866              : 
     867          210 :       NULLIFY (ao_mo_fmstruct)
     868              :       CALL cp_fm_struct_create(ao_mo_fmstruct, nrow_global=nao_red, ncol_global=nmo, &
     869          210 :                                para_env=mo_coeff%matrix_struct%para_env, context=mo_coeff%matrix_struct%context)
     870              : 
     871          210 :       CALL cp_fm_create(u_mo, ao_mo_fmstruct)
     872          210 :       CALL cp_fm_create(u_mo_scaled, ao_mo_fmstruct)
     873              : 
     874          210 :       CALL cp_fm_struct_release(ao_mo_fmstruct)
     875              : 
     876              :       ! U * C
     877          210 :       IF (PRESENT(matrix_u_fm)) THEN
     878          210 :          IF (is_triangular) THEN
     879          124 :             CALL cp_fm_to_fm(mo_coeff, u_mo)
     880              :             CALL cp_fm_triangular_multiply(matrix_u_fm, u_mo, side="L", transpose_tr=.FALSE., &
     881          124 :                                            invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
     882              :          ELSE
     883           86 :             CALL parallel_gemm("N", "N", nao_red, nmo, nao, 1.0_dp, matrix_u_fm, mo_coeff, 0.0_dp, u_mo)
     884              :          END IF
     885              :       ELSE
     886              :          ! assume U is an identity matrix
     887            0 :          CALL cp_fm_to_fm(mo_coeff, u_mo)
     888              :       END IF
     889              : 
     890          210 :       CALL cp_fm_to_fm(u_mo, u_mo_scaled)
     891              : 
     892              :       ! up-shift all unoccupied molecular orbitals by the amount of 'level_shift'
     893              :       ! weight = diag(DELTA) = (0, ... 0, level_shift, ..., level_shift)
     894              :       !             MO index :  1 .. homo   homo+1     ...  nmo
     895          630 :       ALLOCATE (weights(nmo))
     896         1120 :       weights(1:homo) = 0.0_dp
     897         5820 :       weights(homo + 1:nmo) = level_shift
     898              :       ! DELTA * U * C
     899              :       ! DELTA is a diagonal matrix, so simply scale all the columns of (U * C) by weights(:)
     900          210 :       CALL cp_fm_column_scale(u_mo_scaled, weights)
     901          210 :       DEALLOCATE (weights)
     902              : 
     903              :       ! NewKS = U^{-1,T} * KS * U^{-1} + (U * C) * DELTA * (U * C)^T
     904          210 :       CALL parallel_gemm("N", "T", nao_red, nao_red, nmo, 1.0_dp, u_mo, u_mo_scaled, 1.0_dp, matrix_ks_fm)
     905              : 
     906          210 :       CALL cp_fm_release(u_mo_scaled)
     907          210 :       CALL cp_fm_release(u_mo)
     908              : 
     909          210 :       CALL timestop(handle)
     910              : 
     911          420 :    END SUBROUTINE shift_unocc_mos
     912              : 
     913              : END MODULE qs_scf_methods
        

Generated by: LCOV version 2.0-1