LCOV - code coverage report
Current view: top level - src - qs_mom_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 136 173 78.6 %
Date: 2024-03-28 07:31:50 Functions: 4 4 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief methods for deltaSCF calculations
      10             : ! **************************************************************************************************
      11             : MODULE qs_mom_methods
      12             :    USE bibliography,                    ONLY: Barca2018,&
      13             :                                               Gilbert2008,&
      14             :                                               cite_reference
      15             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      16             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      17             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale
      18             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      19             :                                               cp_fm_struct_release,&
      20             :                                               cp_fm_struct_type
      21             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      22             :                                               cp_fm_get_info,&
      23             :                                               cp_fm_maxabsval,&
      24             :                                               cp_fm_to_fm,&
      25             :                                               cp_fm_type,&
      26             :                                               cp_fm_vectorsnorm,&
      27             :                                               cp_fm_vectorssum
      28             :    USE dbcsr_api,                       ONLY: dbcsr_p_type
      29             :    USE input_constants,                 ONLY: momproj_norm,&
      30             :                                               momproj_sum,&
      31             :                                               momtype_imom,&
      32             :                                               momtype_mom
      33             :    USE input_section_types,             ONLY: section_vals_type
      34             :    USE kinds,                           ONLY: dp
      35             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      36             :    USE qs_density_matrices,             ONLY: calculate_density_matrix
      37             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      38             :                                               mo_set_type,&
      39             :                                               set_mo_set
      40             :    USE qs_scf_diagonalization,          ONLY: general_eigenproblem
      41             :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      42             :    USE scf_control_types,               ONLY: scf_control_type
      43             :    USE string_utilities,                ONLY: integer_to_string
      44             :    USE util,                            ONLY: sort,&
      45             :                                               sort_unique
      46             : #include "./base/base_uses.f90"
      47             : 
      48             :    IMPLICIT NONE
      49             : 
      50             :    PRIVATE
      51             : 
      52             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mom_methods'
      53             : 
      54             :    PUBLIC  :: do_mom_guess, do_mom_diag
      55             :    PRIVATE :: mom_is_unique_orbital_indices, mom_reoccupy_orbitals
      56             : 
      57             : CONTAINS
      58             : 
      59             : ! **************************************************************************************************
      60             : !> \brief check that every molecular orbital index appears only once in each
      61             : !>        (de-)occupation list supplied by user. Check that all the indices
      62             : !>        are positive integers and abort if it is not the case.
      63             : !> \param  iarr      list of molecular orbital indices to be checked
      64             : !> \return .true. if all the elements are unique or the list contains
      65             : !>         exactly one 0 element (meaning no excitation)
      66             : !> \par History
      67             : !>      01.2016 created [Sergey Chulkov]
      68             : ! **************************************************************************************************
      69          80 :    FUNCTION mom_is_unique_orbital_indices(iarr) RESULT(is_unique)
      70             :       INTEGER, DIMENSION(:), POINTER                     :: iarr
      71             :       LOGICAL                                            :: is_unique
      72             : 
      73             :       CHARACTER(len=*), PARAMETER :: routineN = 'mom_is_unique_orbital_indices'
      74             : 
      75             :       INTEGER                                            :: handle, norbs
      76             :       INTEGER, DIMENSION(:), POINTER                     :: tmp_iarr
      77             : 
      78          80 :       CALL timeset(routineN, handle)
      79             : 
      80          80 :       CPASSERT(ASSOCIATED(iarr))
      81          80 :       norbs = SIZE(iarr)
      82             : 
      83          80 :       IF (norbs > 0) THEN
      84         240 :          ALLOCATE (tmp_iarr(norbs))
      85             : 
      86         320 :          tmp_iarr(:) = iarr(:)
      87          80 :          CALL sort_unique(tmp_iarr, is_unique)
      88             : 
      89             :          ! Ensure that all orbital indices are positive integers.
      90             :          ! A special value '0' means 'disabled keyword',
      91             :          ! it must appear once to be interpreted in such a way
      92          80 :          IF (tmp_iarr(1) < 0 .OR. (tmp_iarr(1) == 0 .AND. norbs > 1)) &
      93           0 :             CPABORT("MOM: all molecular orbital indices must be positive integer numbers")
      94             : 
      95         160 :          DEALLOCATE (tmp_iarr)
      96             :       END IF
      97             : 
      98             :       is_unique = .TRUE.
      99             : 
     100          80 :       CALL timestop(handle)
     101             : 
     102          80 :    END FUNCTION mom_is_unique_orbital_indices
     103             : 
     104             : ! **************************************************************************************************
     105             : !> \brief swap occupation numbers between molecular orbitals
     106             : !>        from occupation and de-occupation lists
     107             : !> \param mo_set        set of molecular orbitals
     108             : !> \param deocc_orb_set list of de-occupied orbital indices
     109             : !> \param occ_orb_set   list of newly occupied orbital indices
     110             : !> \param spin          spin component of the molecular orbitals;
     111             : !>                      to be used for diagnostic messages
     112             : !> \par History
     113             : !>      01.2016 created [Sergey Chulkov]
     114             : ! **************************************************************************************************
     115          40 :    SUBROUTINE mom_reoccupy_orbitals(mo_set, deocc_orb_set, occ_orb_set, spin)
     116             :       TYPE(mo_set_type), INTENT(INOUT)                   :: mo_set
     117             :       INTEGER, DIMENSION(:), POINTER                     :: deocc_orb_set, occ_orb_set
     118             :       CHARACTER(len=*), INTENT(in)                       :: spin
     119             : 
     120             :       CHARACTER(len=*), PARAMETER :: routineN = 'mom_reoccupy_orbitals'
     121             : 
     122             :       CHARACTER(len=10)                                  :: str_iorb, str_norbs
     123             :       CHARACTER(len=3)                                   :: str_prefix
     124             :       INTEGER                                            :: handle, homo, iorb, lfomo, nao, nmo, &
     125             :                                                             norbs
     126             :       REAL(kind=dp)                                      :: maxocc
     127          40 :       REAL(kind=dp), DIMENSION(:), POINTER               :: occ_nums
     128             : 
     129          40 :       CALL timeset(routineN, handle)
     130             : 
     131             :       ! MOM electron excitation should preserve both the number of electrons and
     132             :       ! multiplicity of the electronic system thus ensuring the following constraint :
     133             :       ! norbs = SIZE(deocc_orb_set) == SIZE(occ_orb_set)
     134          40 :       norbs = SIZE(deocc_orb_set)
     135             : 
     136             :       ! the following assertion should never raise an exception
     137          40 :       CPASSERT(SIZE(deocc_orb_set) == SIZE(occ_orb_set))
     138             : 
     139             :       ! MOM does not follow aufbau principle producing non-uniformly occupied orbitals
     140          40 :       CALL set_mo_set(mo_set=mo_set, uniform_occupation=.FALSE.)
     141             : 
     142          40 :       IF (deocc_orb_set(1) /= 0 .AND. occ_orb_set(1) /= 0) THEN
     143             :          CALL get_mo_set(mo_set=mo_set, maxocc=maxocc, &
     144          20 :                          nao=nao, nmo=nmo, occupation_numbers=occ_nums)
     145             : 
     146          20 :          IF (deocc_orb_set(norbs) > nao .OR. occ_orb_set(norbs) > nao) THEN
     147             :             ! STOP: one of the molecular orbital index exceeds the number of atomic basis functions available
     148           0 :             CALL integer_to_string(nao, str_norbs)
     149             : 
     150           0 :             IF (deocc_orb_set(norbs) >= occ_orb_set(norbs)) THEN
     151           0 :                iorb = deocc_orb_set(norbs)
     152           0 :                str_prefix = 'de-'
     153             :             ELSE
     154           0 :                iorb = occ_orb_set(norbs)
     155           0 :                str_prefix = ''
     156             :             END IF
     157           0 :             CALL integer_to_string(iorb, str_iorb)
     158             : 
     159             :             CALL cp_abort(__LOCATION__, "Unable to "//TRIM(str_prefix)//"occupy "// &
     160             :                           TRIM(spin)//" orbital No. "//TRIM(str_iorb)// &
     161             :                           " since its index exceeds the number of atomic orbital functions available ("// &
     162           0 :                           TRIM(str_norbs)//"). Please consider using a larger basis set.")
     163             :          END IF
     164             : 
     165          20 :          IF (deocc_orb_set(norbs) > nmo .OR. occ_orb_set(norbs) > nmo) THEN
     166             :             ! STOP: one of the molecular orbital index exceeds the number of constructed molecular orbitals
     167           0 :             IF (deocc_orb_set(norbs) >= occ_orb_set(norbs)) THEN
     168           0 :                iorb = deocc_orb_set(norbs)
     169             :             ELSE
     170           0 :                iorb = occ_orb_set(norbs)
     171             :             END IF
     172             : 
     173           0 :             IF (iorb - nmo > 1) THEN
     174           0 :                CALL integer_to_string(iorb - nmo, str_iorb)
     175           0 :                str_prefix = 's'
     176             :             ELSE
     177           0 :                str_iorb = 'an'
     178           0 :                str_prefix = ''
     179             :             END IF
     180             : 
     181           0 :             CALL integer_to_string(nmo, str_norbs)
     182             : 
     183             :             CALL cp_abort(__LOCATION__, "The number of molecular orbitals ("//TRIM(str_norbs)// &
     184             :                           ") is not enough to perform MOM calculation. Please add "// &
     185             :                           TRIM(str_iorb)//" extra orbital"//TRIM(str_prefix)// &
     186           0 :                           " using the ADDED_MOS keyword in the SCF section of your input file.")
     187             :          END IF
     188             : 
     189          40 :          DO iorb = 1, norbs
     190             :             ! swap occupation numbers between two adjoint molecular orbitals
     191          20 :             IF (occ_nums(deocc_orb_set(iorb)) <= 0.0_dp) THEN
     192           0 :                CALL integer_to_string(deocc_orb_set(iorb), str_iorb)
     193             : 
     194             :                CALL cp_abort(__LOCATION__, "The "//TRIM(spin)//" orbital No. "// &
     195           0 :                              TRIM(str_iorb)//" is not occupied thus it cannot be deoccupied.")
     196             :             END IF
     197             : 
     198          20 :             IF (occ_nums(occ_orb_set(iorb)) > 0.0_dp) THEN
     199           0 :                CALL integer_to_string(occ_orb_set(iorb), str_iorb)
     200             : 
     201             :                CALL cp_abort(__LOCATION__, "The "//TRIM(spin)//" orbital No. "// &
     202           0 :                              TRIM(str_iorb)//" is already occupied thus it cannot be reoccupied.")
     203             :             END IF
     204             : 
     205          20 :             occ_nums(occ_orb_set(iorb)) = occ_nums(deocc_orb_set(iorb))
     206          40 :             occ_nums(deocc_orb_set(iorb)) = 0.0_dp
     207             :          END DO
     208             : 
     209             :          ! locate the lowest non-maxocc occupied orbital
     210          78 :          DO lfomo = 1, nmo
     211          78 :             IF (occ_nums(lfomo) /= maxocc) EXIT
     212             :          END DO
     213             : 
     214             :          ! locate the highest occupied orbital
     215          90 :          DO homo = nmo, 1, -1
     216          90 :             IF (occ_nums(homo) > 0.0_dp) EXIT
     217             :          END DO
     218             : 
     219          20 :          CALL set_mo_set(mo_set=mo_set, homo=homo, lfomo=lfomo)
     220             : 
     221          20 :       ELSE IF (deocc_orb_set(1) /= 0 .OR. occ_orb_set(1) /= 0) THEN
     222             :          CALL cp_abort(__LOCATION__, &
     223           0 :                        "Incorrect multiplicity of the MOM reference electronic state")
     224             :       END IF
     225             : 
     226          40 :       CALL timestop(handle)
     227             : 
     228          40 :    END SUBROUTINE mom_reoccupy_orbitals
     229             : 
     230             : ! **************************************************************************************************
     231             : !> \brief initial guess for the maximum overlap method
     232             : !> \param nspins      number of spin components
     233             : !> \param mos         array of molecular orbitals
     234             : !> \param scf_control SCF control variables
     235             : !> \param p_rmpv      density matrix to be computed
     236             : !> \par History
     237             : !>    * 01.2016 created [Sergey Chulkov]
     238             : ! **************************************************************************************************
     239          20 :    SUBROUTINE do_mom_guess(nspins, mos, scf_control, p_rmpv)
     240             :       INTEGER, INTENT(in)                                :: nspins
     241             :       TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT)     :: mos
     242             :       TYPE(scf_control_type), POINTER                    :: scf_control
     243             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: p_rmpv
     244             : 
     245             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'do_mom_guess'
     246             : 
     247             :       CHARACTER(len=10)                                  :: str_iter
     248             :       INTEGER                                            :: handle, ispin, scf_iter
     249             :       LOGICAL                                            :: is_mo
     250             :       REAL(kind=dp)                                      :: maxa
     251             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     252             : 
     253          20 :       CALL timeset(routineN, handle)
     254             : 
     255             :       ! we are about to initialise the maximum overlap method,
     256             :       ! so cite the relevant reference first
     257          20 :       IF (scf_control%diagonalization%mom_type == momtype_mom) THEN
     258          20 :          CALL cite_reference(Gilbert2008)
     259           0 :       ELSE IF (scf_control%diagonalization%mom_type == momtype_imom) THEN
     260           0 :          CALL cite_reference(Barca2018)
     261             :       END IF
     262             : 
     263             :       ! ensure we do not have duplicated orbital indices
     264          20 :       IF (.NOT. &
     265             :           (mom_is_unique_orbital_indices(scf_control%diagonalization%mom_deoccA) .AND. &
     266             :            mom_is_unique_orbital_indices(scf_control%diagonalization%mom_deoccB) .AND. &
     267             :            mom_is_unique_orbital_indices(scf_control%diagonalization%mom_occA) .AND. &
     268             :            mom_is_unique_orbital_indices(scf_control%diagonalization%mom_occB))) &
     269             :          CALL cp_abort(__LOCATION__, &
     270           0 :                        "Duplicate orbital indices were found in the MOM section")
     271             : 
     272             :       ! ignore beta orbitals for spin-unpolarized calculations
     273          20 :       IF (nspins == 1 .AND. (scf_control%diagonalization%mom_deoccB(1) /= 0 &
     274             :                              .OR. scf_control%diagonalization%mom_occB(1) /= 0)) THEN
     275             : 
     276             :          CALL cp_warn(__LOCATION__, "Maximum overlap method will"// &
     277           0 :                       " ignore beta orbitals since neither UKS nor ROKS calculation is performed")
     278             :       END IF
     279             : 
     280             :       ! compute the change in multiplicity and number of electrons
     281             :       IF (SIZE(scf_control%diagonalization%mom_deoccA) /= &
     282          20 :           SIZE(scf_control%diagonalization%mom_occA) .OR. &
     283             :           (nspins > 1 .AND. &
     284             :            SIZE(scf_control%diagonalization%mom_deoccB) /= &
     285             :            SIZE(scf_control%diagonalization%mom_occB))) THEN
     286             : 
     287             :          CALL cp_abort(__LOCATION__, "Incorrect multiplicity of the MOM reference"// &
     288           0 :                        " electronic state or inconsistent number of electrons")
     289             :       END IF
     290             : 
     291          20 :       is_mo = .FALSE.
     292             :       ! by default activate MOM at the second SCF iteration as the
     293             :       ! 'old' molecular orbitals are unavailable from the very beginning
     294          20 :       scf_iter = 2
     295             :       ! check if the molecular orbitals are actually there
     296             :       ! by finding at least one MO coefficient > 0
     297          28 :       DO ispin = 1, nspins
     298          24 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     299          24 :          CALL cp_fm_maxabsval(mo_coeff, maxa)
     300             :          ! is_mo |= maxa > 0.0_dp
     301          28 :          IF (maxa > 0.0_dp) THEN
     302          16 :             is_mo = .TRUE.
     303             :             ! we already have the molecular orbitals (e.g. from a restart file);
     304             :             ! activate MOM immediately if the input keyword START_ITER is not given
     305          16 :             scf_iter = 1
     306          16 :             EXIT
     307             :          END IF
     308             :       END DO
     309             : 
     310             :       ! proceed alpha orbitals
     311          20 :       IF (nspins >= 1) &
     312             :          CALL mom_reoccupy_orbitals(mos(1), &
     313             :                                     scf_control%diagonalization%mom_deoccA, &
     314          20 :                                     scf_control%diagonalization%mom_occA, 'alpha')
     315             : 
     316             :       ! proceed beta orbitals (if any)
     317          20 :       IF (nspins >= 2) &
     318             :          CALL mom_reoccupy_orbitals(mos(2), &
     319             :                                     scf_control%diagonalization%mom_deoccB, &
     320          20 :                                     scf_control%diagonalization%mom_occB, 'beta')
     321             : 
     322             :       ! recompute the density matrix if the molecular orbitals are here;
     323             :       ! otherwise do nothing to prevent zeroing out the density matrix
     324             :       ! obtained from atomic guess
     325          20 :       IF (is_mo) THEN
     326          48 :          DO ispin = 1, nspins
     327          48 :             CALL calculate_density_matrix(mos(ispin), p_rmpv(ispin)%matrix)
     328             :          END DO
     329             :       END IF
     330             : 
     331             :       ! adjust the start SCF iteration number if needed
     332          20 :       IF (scf_control%diagonalization%mom_start < scf_iter) THEN
     333          18 :          IF (scf_control%diagonalization%mom_start > 0) THEN
     334             :             ! inappropriate iteration number has been provided through the input file;
     335             :             ! fix it and issue a warning message
     336           0 :             CALL integer_to_string(scf_iter, str_iter)
     337             :             CALL cp_warn(__LOCATION__, &
     338             :                          "The maximum overlap method will be activated at the SCF iteration No. "// &
     339           0 :                          TRIM(str_iter)//" due to the SCF guess method used.")
     340             :          END IF
     341          18 :          scf_control%diagonalization%mom_start = scf_iter
     342           2 :       ELSE IF (scf_control%diagonalization%mom_start > scf_iter .AND. &
     343             :                (scf_control%diagonalization%mom_occA(1) > 0 .OR. scf_control%diagonalization%mom_occB(1) > 0)) THEN
     344             :          ! the keyword START_ITER has been provided for an excited state calculation, ignore it
     345           2 :          CALL integer_to_string(scf_iter, str_iter)
     346             :          CALL cp_warn(__LOCATION__, &
     347             :                       "The maximum overlap method will be activated at the SCF iteration No. "// &
     348           2 :                       TRIM(str_iter)//" because an excited state calculation has been requested")
     349           2 :          scf_control%diagonalization%mom_start = scf_iter
     350             :       END IF
     351             : 
     352             :       ! MOM is now initialised properly
     353          20 :       scf_control%diagonalization%mom_didguess = .TRUE.
     354             : 
     355          20 :       CALL timestop(handle)
     356             : 
     357          20 :    END SUBROUTINE do_mom_guess
     358             : 
     359             : ! **************************************************************************************************
     360             : !> \brief do an SCF iteration, then compute occupation numbers of the new
     361             : !>  molecular orbitals according to their overlap with the previous ones
     362             : !> \param scf_env     SCF environment information
     363             : !> \param mos         array of molecular orbitals
     364             : !> \param matrix_ks   sparse Kohn-Sham matrix
     365             : !> \param matrix_s    sparse overlap matrix
     366             : !> \param scf_control SCF control variables
     367             : !> \param scf_section SCF input section
     368             : !> \param diis_step   have we done a DIIS step
     369             : !> \par History
     370             : !>    * 07.2014 created [Matt Watkins]
     371             : !>    * 01.2016 release version [Sergey Chulkov]
     372             : !>    * 03.2018 initial maximum overlap method [Sergey Chulkov]
     373             : ! **************************************************************************************************
     374         324 :    SUBROUTINE do_mom_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
     375             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     376             :       TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT)     :: mos
     377             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s
     378             :       TYPE(scf_control_type), POINTER                    :: scf_control
     379             :       TYPE(section_vals_type), POINTER                   :: scf_section
     380             :       LOGICAL, INTENT(INOUT)                             :: diis_step
     381             : 
     382             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'do_mom_diag'
     383             : 
     384             :       INTEGER                                            :: handle, homo, iproj, ispin, lfomo, nao, &
     385             :                                                             nmo, nspins
     386         324 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: inds
     387             :       REAL(kind=dp)                                      :: maxocc
     388         324 :       REAL(kind=dp), DIMENSION(:), POINTER               :: occ_nums, proj, tmp_occ_nums
     389             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     390             :       TYPE(cp_fm_struct_type), POINTER                   :: ao_mo_fmstruct, mo_mo_fmstruct
     391             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff, mo_coeff_ref, overlap, svec
     392             : 
     393         324 :       CALL timeset(routineN, handle)
     394             : 
     395         324 :       IF (.NOT. scf_control%diagonalization%mom_didguess) &
     396             :          CALL cp_abort(__LOCATION__, &
     397           0 :                        "The current implementation of the maximum overlap method is incompatible with the initial SCF guess")
     398             : 
     399             :       ! number of spins == dft_control%nspins
     400         324 :       nspins = SIZE(matrix_ks)
     401             : 
     402             :       ! copy old molecular orbitals
     403         324 :       IF (scf_env%iter_count >= scf_control%diagonalization%mom_start) THEN
     404         318 :          IF (.NOT. ASSOCIATED(scf_env%mom_ref_mo_coeff)) THEN
     405         110 :             ALLOCATE (scf_env%mom_ref_mo_coeff(nspins))
     406          66 :             DO ispin = 1, nspins
     407          44 :                NULLIFY (ao_mo_fmstruct)
     408          44 :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_nums)
     409          44 :                CALL cp_fm_get_info(mo_coeff, matrix_struct=ao_mo_fmstruct)
     410          44 :                CALL cp_fm_create(scf_env%mom_ref_mo_coeff(ispin), ao_mo_fmstruct)
     411             : 
     412             :                ! Initial Maximum Overlap Method: keep initial molecular orbitals
     413          66 :                IF (scf_control%diagonalization%mom_type == momtype_imom) THEN
     414           0 :                   CALL cp_fm_to_fm(mo_coeff, scf_env%mom_ref_mo_coeff(ispin))
     415           0 :                   CALL cp_fm_column_scale(scf_env%mom_ref_mo_coeff(ispin), occ_nums)
     416             :                END IF
     417             :             END DO
     418             :          END IF
     419             : 
     420             :          ! allocate the molecular orbitals overlap matrix
     421         318 :          IF (.NOT. ASSOCIATED(scf_env%mom_overlap)) THEN
     422         100 :             ALLOCATE (scf_env%mom_overlap(nspins))
     423          60 :             DO ispin = 1, nspins
     424          40 :                NULLIFY (blacs_env, mo_mo_fmstruct)
     425          40 :                CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, mo_coeff=mo_coeff)
     426          40 :                CALL cp_fm_get_info(mo_coeff, context=blacs_env)
     427          40 :                CALL cp_fm_struct_create(mo_mo_fmstruct, nrow_global=nmo, ncol_global=nmo, context=blacs_env)
     428          40 :                CALL cp_fm_create(scf_env%mom_overlap(ispin), mo_mo_fmstruct)
     429         100 :                CALL cp_fm_struct_release(mo_mo_fmstruct)
     430             :             END DO
     431             :          END IF
     432             : 
     433             :          ! allocate a matrix to store the product S * mo_coeff
     434         318 :          IF (.NOT. ASSOCIATED(scf_env%mom_s_mo_coeff)) THEN
     435         100 :             ALLOCATE (scf_env%mom_s_mo_coeff(nspins))
     436          60 :             DO ispin = 1, nspins
     437          40 :                NULLIFY (ao_mo_fmstruct)
     438          40 :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     439          40 :                CALL cp_fm_get_info(mo_coeff, matrix_struct=ao_mo_fmstruct)
     440          60 :                CALL cp_fm_create(scf_env%mom_s_mo_coeff(ispin), ao_mo_fmstruct)
     441             :             END DO
     442             :          END IF
     443             : 
     444             :          ! Original Maximum Overlap Method: keep orbitals from the previous SCF iteration
     445         318 :          IF (scf_control%diagonalization%mom_type == momtype_mom) THEN
     446         954 :             DO ispin = 1, nspins
     447         636 :                CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_nums)
     448         636 :                CALL cp_fm_to_fm(mo_coeff, scf_env%mom_ref_mo_coeff(ispin))
     449         954 :                CALL cp_fm_column_scale(scf_env%mom_ref_mo_coeff(ispin), occ_nums)
     450             :             END DO
     451             :          END IF
     452             :       END IF
     453             : 
     454             :       ! solve the eigenproblem
     455         324 :       CALL general_eigenproblem(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
     456             : 
     457         324 :       IF (scf_env%iter_count >= scf_control%diagonalization%mom_start) THEN
     458         954 :          DO ispin = 1, nspins
     459             : 
     460             :             ! TO DO: sparse-matrix variant; check if use_mo_coeff_b is set, and if yes use mo_coeff_b instead
     461             :             CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc, mo_coeff=mo_coeff, &
     462         636 :                             nao=nao, nmo=nmo, occupation_numbers=occ_nums)
     463             : 
     464         636 :             mo_coeff_ref => scf_env%mom_ref_mo_coeff(ispin)
     465         636 :             overlap => scf_env%mom_overlap(ispin)
     466         636 :             svec => scf_env%mom_s_mo_coeff(ispin)
     467             : 
     468             :             ! svec = S * C(new)
     469         636 :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, svec, nmo)
     470             : 
     471             :             ! overlap = C(reference occupied)^T * S * C(new)
     472         636 :             CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff_ref, svec, 0.0_dp, overlap)
     473             : 
     474        1908 :             ALLOCATE (proj(nmo))
     475        1908 :             ALLOCATE (inds(nmo))
     476        1908 :             ALLOCATE (tmp_occ_nums(nmo))
     477             : 
     478             :             ! project the new molecular orbitals into the space of the reference occupied orbitals
     479         636 :             SELECT CASE (scf_control%diagonalization%mom_proj_formula)
     480             :             CASE (momproj_sum)
     481             :                ! proj_j = abs( \sum_i overlap(i, j) )
     482           0 :                CALL cp_fm_vectorssum(overlap, proj)
     483             : 
     484           0 :                DO iproj = 1, nmo
     485           0 :                   proj(iproj) = ABS(proj(iproj))
     486             :                END DO
     487             : 
     488             :             CASE (momproj_norm)
     489             :                ! proj_j = (\sum_i overlap(i, j)**2) ** 0.5
     490         636 :                CALL cp_fm_vectorsnorm(overlap, proj)
     491             : 
     492             :             CASE DEFAULT
     493         636 :                CPABORT("Unimplemented projection formula")
     494             :             END SELECT
     495             : 
     496       11688 :             tmp_occ_nums(:) = occ_nums(:)
     497             :             ! sort occupation numbers in ascending order
     498         636 :             CALL sort(tmp_occ_nums, nmo, inds)
     499             :             ! sort overlap projection in ascending order
     500         636 :             CALL sort(proj, nmo, inds)
     501             : 
     502             :             ! reorder occupation numbers according to overlap projections
     503        5844 :             DO iproj = 1, nmo
     504        5844 :                occ_nums(inds(iproj)) = tmp_occ_nums(iproj)
     505             :             END DO
     506             : 
     507         636 :             DEALLOCATE (tmp_occ_nums)
     508         636 :             DEALLOCATE (inds)
     509         636 :             DEALLOCATE (proj)
     510             : 
     511             :             ! locate the lowest non-fully occupied orbital
     512        2882 :             DO lfomo = 1, nmo
     513        2882 :                IF (occ_nums(lfomo) /= maxocc) EXIT
     514             :             END DO
     515             : 
     516             :             ! locate the highest occupied orbital
     517        2918 :             DO homo = nmo, 1, -1
     518        2918 :                IF (occ_nums(homo) > 0.0_dp) EXIT
     519             :             END DO
     520             : 
     521        1590 :             CALL set_mo_set(mo_set=mos(ispin), homo=homo, lfomo=lfomo)
     522             :          END DO
     523             :       END IF
     524             : 
     525             :       ! recompute density matrix
     526         972 :       DO ispin = 1, nspins
     527         972 :          CALL calculate_density_matrix(mos(ispin), scf_env%p_mix_new(ispin, 1)%matrix)
     528             :       END DO
     529             : 
     530         324 :       CALL timestop(handle)
     531             : 
     532         648 :    END SUBROUTINE do_mom_diag
     533             : 
     534             : END MODULE qs_mom_methods

Generated by: LCOV version 1.15