LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 332 458 72.5 %
Date: 2024-03-28 07:31:50 Functions: 7 7 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             : MODULE qs_tddfpt2_utils
       9             :    USE cell_types,                      ONLY: cell_type
      10             :    USE cp_array_utils,                  ONLY: cp_1d_r_p_type
      11             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      12             :    USE cp_control_types,                ONLY: dft_control_type,&
      13             :                                               tddfpt2_control_type
      14             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      15             :                                               cp_dbcsr_plus_fm_fm_t,&
      16             :                                               cp_dbcsr_sm_fm_multiply,&
      17             :                                               dbcsr_allocate_matrix_set
      18             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_triangular_invert
      19             :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose
      20             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      21             :                                               cp_fm_struct_p_type,&
      22             :                                               cp_fm_struct_release,&
      23             :                                               cp_fm_struct_type
      24             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25             :                                               cp_fm_get_info,&
      26             :                                               cp_fm_release,&
      27             :                                               cp_fm_set_all,&
      28             :                                               cp_fm_to_fm,&
      29             :                                               cp_fm_to_fm_submat,&
      30             :                                               cp_fm_type
      31             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      32             :                                               cp_logger_get_default_io_unit,&
      33             :                                               cp_logger_type
      34             :    USE dbcsr_api,                       ONLY: dbcsr_add,&
      35             :                                               dbcsr_copy,&
      36             :                                               dbcsr_get_info,&
      37             :                                               dbcsr_init_p,&
      38             :                                               dbcsr_p_type,&
      39             :                                               dbcsr_type
      40             :    USE input_constants,                 ONLY: &
      41             :         cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_restore, oe_gllb, oe_lb, oe_none, &
      42             :         oe_saop, oe_shift
      43             :    USE input_section_types,             ONLY: section_vals_create,&
      44             :                                               section_vals_get,&
      45             :                                               section_vals_get_subs_vals,&
      46             :                                               section_vals_release,&
      47             :                                               section_vals_retain,&
      48             :                                               section_vals_set_subs_vals,&
      49             :                                               section_vals_type,&
      50             :                                               section_vals_val_get
      51             :    USE kinds,                           ONLY: dp,&
      52             :                                               int_8
      53             :    USE message_passing,                 ONLY: mp_para_env_type
      54             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      55             :    USE physcon,                         ONLY: evolt
      56             :    USE qs_environment_types,            ONLY: get_qs_env,&
      57             :                                               qs_environment_type
      58             :    USE qs_ks_methods,                   ONLY: qs_ks_build_kohn_sham_matrix
      59             :    USE qs_mo_types,                     ONLY: allocate_mo_set,&
      60             :                                               deallocate_mo_set,&
      61             :                                               get_mo_set,&
      62             :                                               init_mo_set,&
      63             :                                               mo_set_type
      64             :    USE qs_scf_methods,                  ONLY: eigensolver
      65             :    USE qs_scf_post_gpw,                 ONLY: make_lumo_gpw
      66             :    USE qs_scf_types,                    ONLY: ot_method_nr,&
      67             :                                               qs_scf_env_type
      68             :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos
      69             :    USE util,                            ONLY: sort
      70             :    USE xc_pot_saop,                     ONLY: add_saop_pot
      71             :    USE xtb_matrices,                    ONLY: build_xtb_ks_matrix
      72             : #include "./base/base_uses.f90"
      73             : 
      74             :    IMPLICIT NONE
      75             : 
      76             :    PRIVATE
      77             : 
      78             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_utils'
      79             : 
      80             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      81             :    ! number of first derivative components (3: d/dx, d/dy, d/dz)
      82             :    INTEGER, PARAMETER, PRIVATE          :: nderivs = 3
      83             :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
      84             : 
      85             :    PUBLIC :: tddfpt_init_ground_state_mos, tddfpt_release_ground_state_mos
      86             :    PUBLIC :: tddfpt_guess_vectors, tddfpt_init_mos, tddfpt_oecorr
      87             :    PUBLIC :: tddfpt_total_number_of_states
      88             : 
      89             : ! **************************************************************************************************
      90             : 
      91             : CONTAINS
      92             : 
      93             : ! **************************************************************************************************
      94             : !> \brief Prepare MOs for TDDFPT Calculations
      95             : !> \param qs_env  Quickstep environment
      96             : !> \param gs_mos  ...
      97             : !> \param iounit ...
      98             : ! **************************************************************************************************
      99        1054 :    SUBROUTINE tddfpt_init_mos(qs_env, gs_mos, iounit)
     100             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     101             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     102             :          POINTER                                         :: gs_mos
     103             :       INTEGER, INTENT(IN)                                :: iounit
     104             : 
     105             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_init_mos'
     106             : 
     107             :       INTEGER                                            :: handle, ispin, nmo_avail, nmo_occ, &
     108             :                                                             nmo_virt, nspins
     109             :       INTEGER, DIMENSION(2, 2)                           :: moc, mvt
     110             :       LOGICAL                                            :: print_virtuals_newtonx
     111        1054 :       REAL(kind=dp), DIMENSION(:), POINTER               :: evals_virt_spin
     112             :       TYPE(cell_type), POINTER                           :: cell
     113        1054 :       TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER        :: evals_virt
     114             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     115             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
     116        1054 :          TARGET                                          :: mos_virt
     117             :       TYPE(cp_fm_type), POINTER                          :: mos_virt_spin
     118        1054 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s
     119             :       TYPE(dft_control_type), POINTER                    :: dft_control
     120        1054 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     121             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     122             :       TYPE(section_vals_type), POINTER                   :: print_section
     123             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     124             : 
     125        1054 :       CALL timeset(routineN, handle)
     126             : 
     127             :       CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
     128        1054 :                       matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
     129        1054 :       tddfpt_control => dft_control%tddfpt2_control
     130             : 
     131        1054 :       CPASSERT(.NOT. ASSOCIATED(gs_mos))
     132             :       ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
     133        1054 :       nspins = dft_control%nspins
     134        4340 :       ALLOCATE (gs_mos(nspins))
     135             : 
     136             :       ! check if virtuals should be constructed for NAMD interface with NEWTONX
     137        1054 :       print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
     138        1054 :       CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals_newtonx)
     139             : 
     140             :       ! when the number of unoccupied orbitals is limited and OT has been used
     141             :       ! for the ground-state DFT calculation,
     142             :       ! compute the missing unoccupied orbitals using OT as well.
     143        1054 :       NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
     144        1054 :       IF (ASSOCIATED(scf_env)) THEN
     145        1054 :          IF ((scf_env%method == ot_method_nr .AND. tddfpt_control%nlumo > 0) .OR. &
     146             :              (scf_env%method == ot_method_nr .AND. print_virtuals_newtonx)) THEN
     147             :             ! As OT with ADDED_MOS/=0 is currently not implemented, the following block is equivalent to:
     148             :             ! nmo_virt = tddfpt_control%nlumo
     149             :             ! number of already computed unoccupied orbitals (added_mos) .
     150           2 :             nmo_virt = HUGE(0)
     151           4 :             DO ispin = 1, nspins
     152           2 :                CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
     153           4 :                nmo_virt = MIN(nmo_virt, nmo_avail - nmo_occ)
     154             :             END DO
     155             :             ! number of unoccupied orbitals to compute
     156           2 :             nmo_virt = tddfpt_control%nlumo - nmo_virt
     157           2 :             IF (.NOT. print_virtuals_newtonx) THEN
     158           0 :                IF (nmo_virt > 0) THEN
     159           0 :                   ALLOCATE (evals_virt(nspins), mos_virt(nspins))
     160             :                   ! the number of actually computed unoccupied orbitals will be stored as 'nmo_avail'
     161           0 :                   CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
     162             :                END IF
     163             :             END IF
     164             :          END IF
     165             :       END IF
     166             : 
     167        2232 :       DO ispin = 1, nspins
     168        1178 :          IF (ASSOCIATED(evals_virt)) THEN
     169           0 :             evals_virt_spin => evals_virt(ispin)%array
     170             :          ELSE
     171        1178 :             NULLIFY (evals_virt_spin)
     172             :          END IF
     173        1178 :          IF (ALLOCATED(mos_virt)) THEN
     174           0 :             mos_virt_spin => mos_virt(ispin)
     175             :          ELSE
     176        1178 :             NULLIFY (mos_virt_spin)
     177             :          END IF
     178             :          CALL tddfpt_init_ground_state_mos(gs_mos=gs_mos(ispin), mo_set=mos(ispin), &
     179             :                                            nlumo=tddfpt_control%nlumo, &
     180             :                                            blacs_env=blacs_env, cholesky_method=cholesky_restore, &
     181             :                                            matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
     182             :                                            mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, &
     183        2232 :                                            qs_env=qs_env)
     184             :       END DO
     185             : 
     186        1054 :       moc = 0
     187        1054 :       mvt = 0
     188        2232 :       DO ispin = 1, nspins
     189        1178 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
     190        2232 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
     191             :       END DO
     192        1054 :       IF (iounit > 0) THEN
     193         527 :          WRITE (iounit, "(T2,A,T36,A)") "TDDFPT| Molecular Orbitals:", &
     194        1054 :             " Spin       AOs       Occ      Virt     Total"
     195        1116 :          DO ispin = 1, nspins
     196         589 :             WRITE (iounit, "(T2,A,T37,I4,4I10)") "TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
     197        1705 :                mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
     198             :          END DO
     199             :       END IF
     200             : 
     201        1054 :       IF (ASSOCIATED(evals_virt)) THEN
     202           0 :          DO ispin = 1, SIZE(evals_virt)
     203           0 :             IF (ASSOCIATED(evals_virt(ispin)%array)) DEALLOCATE (evals_virt(ispin)%array)
     204             :          END DO
     205           0 :          DEALLOCATE (evals_virt)
     206             :       END IF
     207             : 
     208        1054 :       CALL cp_fm_release(mos_virt)
     209             : 
     210        1054 :       CALL timestop(handle)
     211             : 
     212        3162 :    END SUBROUTINE tddfpt_init_mos
     213             : 
     214             : ! **************************************************************************************************
     215             : !> \brief Generate all virtual molecular orbitals for a given spin by diagonalising
     216             : !>        the corresponding Kohn-Sham matrix.
     217             : !> \param gs_mos           structure to store occupied and virtual molecular orbitals
     218             : !>                         (allocated and initialised on exit)
     219             : !> \param mo_set           ground state molecular orbitals for a given spin
     220             : !> \param nlumo            number of unoccupied states to consider (-1 means all states)
     221             : !> \param blacs_env        BLACS parallel environment
     222             : !> \param cholesky_method  Cholesky method to compute the inverse overlap matrix
     223             : !> \param matrix_ks        Kohn-Sham matrix for a given spin
     224             : !> \param matrix_s         overlap matrix
     225             : !> \param mos_virt         precomputed (OT) expansion coefficients of virtual molecular orbitals
     226             : !>                         (in addition to the ADDED_MOS, if present). NULL when no OT is in use.
     227             : !> \param evals_virt       orbital energies of precomputed (OT) virtual molecular orbitals.
     228             : !>                         NULL when no OT is in use.
     229             : !> \param qs_env ...
     230             : !> \par History
     231             : !>    * 05.2016 created as tddfpt_lumos() [Sergey Chulkov]
     232             : !>    * 06.2016 renamed, altered prototype [Sergey Chulkov]
     233             : !>    * 04.2019 limit the number of unoccupied states, orbital energy correction [Sergey Chulkov]
     234             : ! **************************************************************************************************
     235        1178 :    SUBROUTINE tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, &
     236             :                                            mos_virt, evals_virt, qs_env)
     237             :       TYPE(tddfpt_ground_state_mos)                      :: gs_mos
     238             :       TYPE(mo_set_type), INTENT(IN)                      :: mo_set
     239             :       INTEGER, INTENT(in)                                :: nlumo
     240             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     241             :       INTEGER, INTENT(in)                                :: cholesky_method
     242             :       TYPE(dbcsr_type), POINTER                          :: matrix_ks, matrix_s
     243             :       TYPE(cp_fm_type), INTENT(IN), POINTER              :: mos_virt
     244             :       REAL(kind=dp), DIMENSION(:), POINTER               :: evals_virt
     245             :       TYPE(qs_environment_type), INTENT(in), POINTER     :: qs_env
     246             : 
     247             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_init_ground_state_mos'
     248             :       REAL(kind=dp), PARAMETER                           :: eps_dp = EPSILON(0.0_dp)
     249             : 
     250             :       INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, &
     251             :          irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int
     252        1178 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: minrow_neg_array, minrow_pos_array, &
     253        1178 :                                                             sum_sign_array
     254        1178 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     255             :       LOGICAL                                            :: do_eigen, print_phases
     256             :       REAL(kind=dp)                                      :: element, maxocc
     257             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
     258        1178 :          POINTER                                         :: my_block
     259        1178 :       REAL(kind=dp), DIMENSION(:), POINTER               :: mo_evals_extended, mo_occ_extended, &
     260        1178 :                                                             mo_occ_scf
     261             :       TYPE(cp_fm_struct_type), POINTER                   :: ao_ao_fm_struct, ao_mo_occ_fm_struct, &
     262             :                                                             ao_mo_virt_fm_struct, wfn_fm_struct
     263             :       TYPE(cp_fm_type)                                   :: matrix_ks_fm, ortho_fm, work_fm, &
     264             :                                                             work_fm_virt
     265             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff_extended
     266             :       TYPE(cp_logger_type), POINTER                      :: logger
     267             :       TYPE(mo_set_type), POINTER                         :: mos_extended
     268             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     269             :       TYPE(section_vals_type), POINTER                   :: print_section
     270             : 
     271        1178 :       CALL timeset(routineN, handle)
     272             : 
     273        1178 :       NULLIFY (logger)
     274        1178 :       logger => cp_get_default_logger()
     275        1178 :       iounit = cp_logger_get_default_io_unit(logger)
     276             : 
     277        1178 :       CALL blacs_env%get(para_env=para_env)
     278             : 
     279             :       CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, &
     280        1178 :                       nelectron=nelectrons, occupation_numbers=mo_occ_scf)
     281             : 
     282        1178 :       print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
     283        1178 :       CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
     284             : 
     285        1178 :       nmo_virt = nao - nmo_occ
     286        1178 :       IF (nlumo >= 0) &
     287           0 :          nmo_virt = MIN(nmo_virt, nlumo)
     288             : 
     289        1178 :       IF (nmo_virt <= 0) &
     290             :          CALL cp_abort(__LOCATION__, &
     291           0 :                        'At least one unoccupied molecular orbital is required to calculate excited states.')
     292             : 
     293        1178 :       do_eigen = .FALSE.
     294             :       ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
     295        1178 :       IF (ASSOCIATED(evals_virt)) THEN
     296           0 :          CPASSERT(ASSOCIATED(mos_virt))
     297           0 :          IF (nmo_virt > nmo_scf - nmo_occ + SIZE(evals_virt)) do_eigen = .TRUE.
     298             :       ELSE
     299        1178 :          IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .TRUE.
     300             :       END IF
     301             : 
     302             :       ! ++ allocate storage space for gs_mos
     303        1178 :       NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
     304        1178 :       CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
     305        1178 :       CALL cp_fm_struct_create(ao_mo_virt_fm_struct, nrow_global=nao, ncol_global=nmo_virt, context=blacs_env)
     306             : 
     307        1178 :       NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
     308        1178 :       ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
     309        1178 :       CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct)
     310        1178 :       CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
     311             : 
     312        3534 :       ALLOCATE (gs_mos%evals_occ(nmo_occ))
     313        3534 :       ALLOCATE (gs_mos%evals_virt(nmo_virt))
     314        3534 :       ALLOCATE (gs_mos%phases_occ(nmo_occ))
     315        3534 :       ALLOCATE (gs_mos%phases_virt(nmo_virt))
     316             : 
     317             :       ! ++ nullify pointers
     318        1178 :       NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
     319        1178 :       NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
     320        1178 :       CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
     321             : 
     322        1178 :       IF (do_eigen) THEN
     323             :          ! ++ set of molecular orbitals
     324        1178 :          CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
     325        1178 :          ALLOCATE (mos_extended)
     326             :          CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
     327        1178 :                               REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp)
     328        1178 :          CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name="mos-extended")
     329        1178 :          CALL cp_fm_struct_release(wfn_fm_struct)
     330             :          CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
     331        1178 :                          eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
     332             : 
     333             :          ! use the explicit loop in order to avoid temporary arrays.
     334             :          !
     335             :          ! The assignment statement : mo_occ_extended(1:nmo_scf) = mo_occ_scf(1:nmo_scf)
     336             :          ! implies temporary arrays as a compiler does not know in advance that the pointers
     337             :          ! on both sides of the statement point to non-overlapped memory regions
     338        7448 :          DO imo = 1, nmo_scf
     339        7448 :             mo_occ_extended(imo) = mo_occ_scf(imo)
     340             :          END DO
     341       20430 :          mo_occ_extended(nmo_scf + 1:) = 0.0_dp
     342             : 
     343             :          ! ++ allocate temporary matrices
     344        1178 :          CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct)
     345        1178 :          CALL cp_fm_create(ortho_fm, ao_ao_fm_struct)
     346        1178 :          CALL cp_fm_create(work_fm, ao_ao_fm_struct)
     347        1178 :          CALL cp_fm_struct_release(ao_ao_fm_struct)
     348             : 
     349             :          ! some stuff from the subroutine general_eigenproblem()
     350        1178 :          CALL copy_dbcsr_to_fm(matrix_s, ortho_fm)
     351        1178 :          CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm)
     352             : 
     353        1178 :          IF (cholesky_method == cholesky_dbcsr) THEN
     354           0 :             CPABORT('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
     355        1178 :          ELSE IF (cholesky_method == cholesky_off) THEN
     356           0 :             CPABORT('CHOLESKY OFF is not implemented in TDDFT.')
     357             :          ELSE
     358        1178 :             CALL cp_fm_cholesky_decompose(ortho_fm)
     359        1178 :             IF (cholesky_method == cholesky_inverse) THEN
     360           0 :                CALL cp_fm_triangular_invert(ortho_fm)
     361             :             END IF
     362             : 
     363             :             ! need to store 'cholesky_method' in a temporary variable, as the subroutine eigensolver()
     364             :             ! will update this variable
     365        1178 :             cholesky_method_inout = cholesky_method
     366             :             CALL eigensolver(matrix_ks_fm=matrix_ks_fm, mo_set=mos_extended, ortho=ortho_fm, &
     367             :                              work=work_fm, cholesky_method=cholesky_method_inout, &
     368        1178 :                              do_level_shift=.FALSE., level_shift=0.0_dp, use_jacobi=.FALSE.)
     369             :          END IF
     370             : 
     371             :          ! -- clean up needless matrices
     372        1178 :          CALL cp_fm_release(work_fm)
     373        1178 :          CALL cp_fm_release(ortho_fm)
     374        1178 :          CALL cp_fm_release(matrix_ks_fm)
     375             :       ELSE
     376             :          CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
     377           0 :                          eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
     378             :       END IF
     379             : 
     380             :       ! compute the phase of molecular orbitals;
     381             :       ! matrix work_fm holds occupied molecular orbital coefficients distributed among all the processors
     382             :       !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
     383        1178 :       CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
     384        1178 :       CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
     385             : 
     386        1178 :       CALL cp_fm_to_fm(mo_coeff_extended, work_fm, ncol=nmo_occ, source_start=1, target_start=1)
     387             :       CALL cp_fm_get_info(work_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
     388        1178 :                           row_indices=row_indices, col_indices=col_indices, local_data=my_block)
     389             : 
     390        8246 :       ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
     391        7448 :       minrow_neg_array(:) = nao
     392        7448 :       minrow_pos_array(:) = nao
     393        7448 :       sum_sign_array(:) = 0
     394        7448 :       DO icol_local = 1, ncol_local
     395        6270 :          icol_global = col_indices(icol_local)
     396             : 
     397      222770 :          DO irow_local = 1, nrow_local
     398      215322 :             element = my_block(irow_local, icol_local)
     399             : 
     400      215322 :             sign_int = 0
     401      215322 :             IF (element >= eps_dp) THEN
     402             :                sign_int = 1
     403      108805 :             ELSE IF (element <= -eps_dp) THEN
     404      108256 :                sign_int = -1
     405             :             END IF
     406             : 
     407      215322 :             sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
     408             : 
     409      215322 :             irow_global = row_indices(irow_local)
     410      221592 :             IF (sign_int > 0) THEN
     411      106517 :                IF (minrow_pos_array(icol_global) > irow_global) &
     412        6162 :                   minrow_pos_array(icol_global) = irow_global
     413      108805 :             ELSE IF (sign_int < 0) THEN
     414      108256 :                IF (minrow_neg_array(icol_global) > irow_global) &
     415        6018 :                   minrow_neg_array(icol_global) = irow_global
     416             :             END IF
     417             :          END DO
     418             :       END DO
     419             : 
     420        1178 :       CALL para_env%sum(sum_sign_array)
     421        1178 :       CALL para_env%min(minrow_neg_array)
     422        1178 :       CALL para_env%min(minrow_pos_array)
     423             : 
     424        7448 :       DO icol_local = 1, nmo_occ
     425        7448 :          IF (sum_sign_array(icol_local) > 0) THEN
     426             :             ! most of the expansion coefficients are positive => MO's phase = +1
     427        2770 :             gs_mos%phases_occ(icol_local) = 1.0_dp
     428        3500 :          ELSE IF (sum_sign_array(icol_local) < 0) THEN
     429             :             ! most of the expansion coefficients are negative => MO's phase = -1
     430        3254 :             gs_mos%phases_occ(icol_local) = -1.0_dp
     431             :          ELSE
     432             :             ! equal number of positive and negative expansion coefficients
     433         246 :             IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
     434             :                ! the first positive expansion coefficient has a lower index then
     435             :                ! the first negative expansion coefficient; MO's phase = +1
     436         128 :                gs_mos%phases_occ(icol_local) = 1.0_dp
     437             :             ELSE
     438             :                ! MO's phase = -1
     439         118 :                gs_mos%phases_occ(icol_local) = -1.0_dp
     440             :             END IF
     441             :          END IF
     442             :       END DO
     443             : 
     444        1178 :       DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
     445             : 
     446             :       ! return the requested occupied and virtual molecular orbitals and corresponding orbital energies
     447        1178 :       CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
     448        7448 :       gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
     449             : 
     450        1178 :       IF (ASSOCIATED(evals_virt) .AND. (.NOT. do_eigen) .AND. nmo_virt > nmo_scf - nmo_occ) THEN
     451             :          CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_scf - nmo_occ, &
     452           0 :                           source_start=nmo_occ + 1, target_start=1)
     453             :          CALL cp_fm_to_fm(mos_virt, gs_mos%mos_virt, ncol=nmo_virt - (nmo_scf - nmo_occ), &
     454           0 :                           source_start=1, target_start=nmo_scf - nmo_occ + 1)
     455           0 :          gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
     456           0 :          gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
     457             :       ELSE
     458        1178 :          CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
     459       20430 :          gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
     460             :       END IF
     461             : 
     462        1178 :       IF (print_phases) THEN
     463             :          ! compute the phase of molecular orbitals;
     464             :          ! matrix work_fm holds virtual molecular orbital coefficients distributed among all the processors
     465             :          !CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
     466           0 :          CALL cp_fm_create(work_fm_virt, ao_mo_virt_fm_struct)
     467             : 
     468           0 :          CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1)
     469             :          CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, &
     470           0 :                              row_indices=row_indices, col_indices=col_indices, local_data=my_block)
     471             : 
     472           0 :          ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
     473           0 :          minrow_neg_array(:) = nao
     474           0 :          minrow_pos_array(:) = nao
     475           0 :          sum_sign_array(:) = 0
     476           0 :          DO icol_local = 1, ncol_local
     477           0 :             icol_global = col_indices(icol_local)
     478             : 
     479           0 :             DO irow_local = 1, nrow_local
     480           0 :                element = my_block(irow_local, icol_local)
     481             : 
     482           0 :                sign_int = 0
     483           0 :                IF (element >= eps_dp) THEN
     484             :                   sign_int = 1
     485           0 :                ELSE IF (element <= -eps_dp) THEN
     486           0 :                   sign_int = -1
     487             :                END IF
     488             : 
     489           0 :                sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
     490             : 
     491           0 :                irow_global = row_indices(irow_local)
     492           0 :                IF (sign_int > 0) THEN
     493           0 :                   IF (minrow_pos_array(icol_global) > irow_global) &
     494           0 :                      minrow_pos_array(icol_global) = irow_global
     495           0 :                ELSE IF (sign_int < 0) THEN
     496           0 :                   IF (minrow_neg_array(icol_global) > irow_global) &
     497           0 :                      minrow_neg_array(icol_global) = irow_global
     498             :                END IF
     499             :             END DO
     500             :          END DO
     501             : 
     502           0 :          CALL para_env%sum(sum_sign_array)
     503           0 :          CALL para_env%min(minrow_neg_array)
     504           0 :          CALL para_env%min(minrow_pos_array)
     505           0 :          DO icol_local = 1, nmo_virt
     506           0 :             IF (sum_sign_array(icol_local) > 0) THEN
     507             :                ! most of the expansion coefficients are positive => MO's phase = +1
     508           0 :                gs_mos%phases_virt(icol_local) = 1.0_dp
     509           0 :             ELSE IF (sum_sign_array(icol_local) < 0) THEN
     510             :                ! most of the expansion coefficients are negative => MO's phase = -1
     511           0 :                gs_mos%phases_virt(icol_local) = -1.0_dp
     512             :             ELSE
     513             :                ! equal number of positive and negative expansion coefficients
     514           0 :                IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local)) THEN
     515             :                   ! the first positive expansion coefficient has a lower index then
     516             :                   ! the first negative expansion coefficient; MO's phase = +1
     517           0 :                   gs_mos%phases_virt(icol_local) = 1.0_dp
     518             :                ELSE
     519             :                   ! MO's phase = -1
     520           0 :                   gs_mos%phases_virt(icol_local) = -1.0_dp
     521             :                END IF
     522             :             END IF
     523             :          END DO
     524           0 :          DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
     525           0 :          CALL cp_fm_release(work_fm_virt)
     526             :       END IF !print_phases
     527        1178 :       CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
     528             : 
     529        1178 :       CALL cp_fm_release(work_fm)
     530             : 
     531        1178 :       IF (do_eigen) THEN
     532        1178 :          CALL deallocate_mo_set(mos_extended)
     533        1178 :          DEALLOCATE (mos_extended)
     534             :       END IF
     535             : 
     536        1178 :       CALL timestop(handle)
     537             : 
     538        5890 :    END SUBROUTINE tddfpt_init_ground_state_mos
     539             : 
     540             : ! **************************************************************************************************
     541             : !> \brief Release molecular orbitals.
     542             : !> \param gs_mos  structure that holds occupied and virtual molecular orbitals
     543             : !> \par History
     544             : !>    * 06.2016 created [Sergey Chulkov]
     545             : ! **************************************************************************************************
     546        1178 :    SUBROUTINE tddfpt_release_ground_state_mos(gs_mos)
     547             :       TYPE(tddfpt_ground_state_mos), INTENT(inout)       :: gs_mos
     548             : 
     549             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_release_ground_state_mos'
     550             : 
     551             :       INTEGER                                            :: handle
     552             : 
     553        1178 :       CALL timeset(routineN, handle)
     554             : 
     555        1178 :       IF (ALLOCATED(gs_mos%phases_occ)) &
     556        1178 :          DEALLOCATE (gs_mos%phases_occ)
     557             : 
     558        1178 :       IF (ALLOCATED(gs_mos%evals_virt)) &
     559        1178 :          DEALLOCATE (gs_mos%evals_virt)
     560             : 
     561        1178 :       IF (ALLOCATED(gs_mos%evals_occ)) &
     562        1178 :          DEALLOCATE (gs_mos%evals_occ)
     563             : 
     564        1178 :       IF (ALLOCATED(gs_mos%phases_virt)) &
     565        1178 :          DEALLOCATE (gs_mos%phases_virt)
     566             : 
     567        1178 :       IF (ASSOCIATED(gs_mos%evals_occ_matrix)) THEN
     568          36 :          CALL cp_fm_release(gs_mos%evals_occ_matrix)
     569          36 :          DEALLOCATE (gs_mos%evals_occ_matrix)
     570             :       END IF
     571             : 
     572        1178 :       IF (ASSOCIATED(gs_mos%mos_virt)) THEN
     573        1178 :          CALL cp_fm_release(gs_mos%mos_virt)
     574        1178 :          DEALLOCATE (gs_mos%mos_virt)
     575             :       END IF
     576             : 
     577        1178 :       IF (ASSOCIATED(gs_mos%mos_occ)) THEN
     578        1178 :          CALL cp_fm_release(gs_mos%mos_occ)
     579        1178 :          DEALLOCATE (gs_mos%mos_occ)
     580             :       END IF
     581             : 
     582        1178 :       CALL timestop(handle)
     583        1178 :    END SUBROUTINE tddfpt_release_ground_state_mos
     584             : 
     585             : ! **************************************************************************************************
     586             : !> \brief Callculate orbital corrected KS matrix for TDDFPT
     587             : !> \param qs_env  Quickstep environment
     588             : !> \param gs_mos ...
     589             : !> \param matrix_ks_oep ...
     590             : ! **************************************************************************************************
     591        1054 :    SUBROUTINE tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
     592             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     593             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     594             :          POINTER                                         :: gs_mos
     595             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks_oep
     596             : 
     597             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_oecorr'
     598             : 
     599             :       INTEGER                                            :: handle, iounit, ispin, nao, nmo_occ, &
     600             :                                                             nspins
     601             :       LOGICAL                                            :: do_hfx
     602             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     603             :       TYPE(cp_fm_struct_type), POINTER                   :: ao_mo_occ_fm_struct, &
     604             :                                                             mo_occ_mo_occ_fm_struct
     605             :       TYPE(cp_fm_type)                                   :: work_fm
     606             :       TYPE(cp_logger_type), POINTER                      :: logger
     607        1054 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     608             :       TYPE(dft_control_type), POINTER                    :: dft_control
     609             :       TYPE(section_vals_type), POINTER                   :: hfx_section, xc_fun_empty, &
     610             :                                                             xc_fun_original
     611             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     612             : 
     613        1054 :       CALL timeset(routineN, handle)
     614             : 
     615        1054 :       NULLIFY (logger)
     616        1054 :       logger => cp_get_default_logger()
     617        1054 :       iounit = cp_logger_get_default_io_unit(logger)
     618             : 
     619        1054 :       CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
     620        1054 :       tddfpt_control => dft_control%tddfpt2_control
     621             : 
     622             :       ! obtain corrected KS-matrix
     623             :       ! We should 'save' the energy values?
     624        1054 :       nspins = SIZE(gs_mos)
     625        1054 :       NULLIFY (matrix_ks_oep)
     626        1054 :       IF (tddfpt_control%oe_corr /= oe_none) THEN
     627          30 :          IF (iounit > 0) THEN
     628          15 :             WRITE (iounit, "(1X,A)") "", &
     629          15 :                "-------------------------------------------------------------------------------", &
     630          15 :                "-                    Orbital Eigenvalue Correction Started                    -", &
     631          30 :                "-------------------------------------------------------------------------------"
     632             :          END IF
     633             : 
     634             :          CALL cp_warn(__LOCATION__, &
     635             :                       "Orbital energy correction potential is an experimental feature. "// &
     636          30 :                       "Use it with extreme care")
     637             : 
     638          30 :          hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
     639          30 :          CALL section_vals_get(hfx_section, explicit=do_hfx)
     640          30 :          IF (do_hfx) THEN
     641             :             CALL cp_abort(__LOCATION__, &
     642             :                           "Implementation of orbital energy correction XC-potentials is "// &
     643           0 :                           "currently incompatible with exact-exchange functionals")
     644             :          END IF
     645             : 
     646          30 :          CALL dbcsr_allocate_matrix_set(matrix_ks_oep, nspins)
     647          66 :          DO ispin = 1, nspins
     648          36 :             CALL dbcsr_init_p(matrix_ks_oep(ispin)%matrix)
     649          66 :             CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
     650             :          END DO
     651             : 
     652             :          ! KS-matrix without XC-terms
     653          30 :          xc_fun_original => section_vals_get_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL")
     654          30 :          CALL section_vals_retain(xc_fun_original)
     655          30 :          NULLIFY (xc_fun_empty)
     656          30 :          CALL section_vals_create(xc_fun_empty, xc_fun_original%section)
     657          30 :          CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_empty)
     658          30 :          CALL section_vals_release(xc_fun_empty)
     659             : 
     660          30 :          IF (dft_control%qs_control%semi_empirical) THEN
     661           0 :             CPABORT("TDDFPT with SE not possible")
     662          30 :          ELSEIF (dft_control%qs_control%dftb) THEN
     663           0 :             CPABORT("TDDFPT with DFTB not possible")
     664          30 :          ELSEIF (dft_control%qs_control%xtb) THEN
     665             :             CALL build_xtb_ks_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
     666          16 :                                      ext_ks_matrix=matrix_ks_oep)
     667             :          ELSE
     668             :             CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
     669          14 :                                               ext_ks_matrix=matrix_ks_oep)
     670             :          END IF
     671             : 
     672             :          IF (tddfpt_control%oe_corr == oe_saop .OR. &
     673          30 :              tddfpt_control%oe_corr == oe_lb .OR. &
     674             :              tddfpt_control%oe_corr == oe_gllb) THEN
     675          14 :             IF (iounit > 0) THEN
     676           7 :                WRITE (iounit, "(T2,A)") " Orbital energy correction of SAOP type "
     677             :             END IF
     678          14 :             CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
     679          16 :          ELSE IF (tddfpt_control%oe_corr == oe_shift) THEN
     680          16 :             IF (iounit > 0) THEN
     681             :                WRITE (iounit, "(T2,A,T71,F10.3)") &
     682           8 :                   " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*evolt
     683             :                WRITE (iounit, "(T2,A,T71,F10.3)") &
     684           8 :                   " Open Shell Orbital Eigenvalue Shift [eV] ", tddfpt_control%eos_shift*evolt
     685             :             END IF
     686             :             CALL ev_shift_operator(qs_env, gs_mos, matrix_ks_oep, &
     687          16 :                                    tddfpt_control%ev_shift, tddfpt_control%eos_shift)
     688             :          ELSE
     689             :             CALL cp_abort(__LOCATION__, &
     690           0 :                           "Unimplemented orbital energy correction potential")
     691             :          END IF
     692          30 :          CALL section_vals_set_subs_vals(qs_env%input, "DFT%XC%XC_FUNCTIONAL", xc_fun_original)
     693          30 :          CALL section_vals_release(xc_fun_original)
     694             : 
     695             :          ! compute 'evals_occ_matrix'
     696          30 :          CALL dbcsr_get_info(matrix_ks(1)%matrix, nfullrows_total=nao)
     697          30 :          NULLIFY (mo_occ_mo_occ_fm_struct)
     698          66 :          DO ispin = 1, nspins
     699          36 :             nmo_occ = SIZE(gs_mos(ispin)%evals_occ)
     700             :             CALL cp_fm_struct_create(mo_occ_mo_occ_fm_struct, nrow_global=nmo_occ, ncol_global=nmo_occ, &
     701          36 :                                      context=blacs_env)
     702          36 :             ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
     703          36 :             CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
     704          36 :             CALL cp_fm_struct_release(mo_occ_mo_occ_fm_struct)
     705             :             ! work_fm is a temporary [nao x nmo_occ] matrix
     706             :             CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, &
     707          36 :                                      context=blacs_env)
     708          36 :             CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
     709          36 :             CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
     710             :             CALL cp_dbcsr_sm_fm_multiply(matrix_ks_oep(ispin)%matrix, gs_mos(ispin)%mos_occ, &
     711          36 :                                          work_fm, ncol=nmo_occ, alpha=1.0_dp, beta=0.0_dp)
     712             :             CALL parallel_gemm('T', 'N', nmo_occ, nmo_occ, nao, 1.0_dp, gs_mos(ispin)%mos_occ, work_fm, &
     713          36 :                                0.0_dp, gs_mos(ispin)%evals_occ_matrix)
     714          66 :             CALL cp_fm_release(work_fm)
     715             :          END DO
     716          30 :          IF (iounit > 0) THEN
     717             :             WRITE (iounit, "(1X,A)") &
     718          15 :                "-------------------------------------------------------------------------------"
     719             :          END IF
     720             : 
     721             :       END IF
     722             : 
     723        1054 :       CALL timestop(handle)
     724             : 
     725        1054 :    END SUBROUTINE tddfpt_oecorr
     726             : 
     727             : ! **************************************************************************************************
     728             : !> \brief Compute the number of possible singly excited states (occ -> virt)
     729             : !> \param gs_mos          occupied and virtual molecular orbitals optimised for the ground state
     730             : !> \return the number of possible single excitations
     731             : !> \par History
     732             : !>    * 01.2017 created [Sergey Chulkov]
     733             : ! **************************************************************************************************
     734        2653 :    PURE FUNCTION tddfpt_total_number_of_states(gs_mos) RESULT(nstates_total)
     735             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     736             :          INTENT(in)                                      :: gs_mos
     737             :       INTEGER(kind=int_8)                                :: nstates_total
     738             : 
     739             :       INTEGER                                            :: ispin, nspins
     740             : 
     741        2653 :       nstates_total = 0
     742        2653 :       nspins = SIZE(gs_mos)
     743             : 
     744        5614 :       DO ispin = 1, nspins
     745             :          nstates_total = nstates_total + &
     746             :                          SIZE(gs_mos(ispin)%evals_occ, kind=int_8)* &
     747        5614 :                          SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
     748             :       END DO
     749        2653 :    END FUNCTION tddfpt_total_number_of_states
     750             : 
     751             : ! **************************************************************************************************
     752             : !> \brief Create a shift operator on virtual/open shell space
     753             : !>        Shift operator = Edelta*Q  Q: projector on virtual space (1-PS)
     754             : !>                                      projector on open shell space PosS
     755             : !> \param qs_env the qs_env that is perturbed by this p_env
     756             : !> \param gs_mos  ...
     757             : !> \param matrix_ks ...
     758             : !> \param ev_shift ...
     759             : !> \param eos_shift ...
     760             : !> \par History
     761             : !>      02.04.2019 adapted for TDDFT use from p_env (JGH)
     762             : !> \author JGH
     763             : ! **************************************************************************************************
     764          16 :    SUBROUTINE ev_shift_operator(qs_env, gs_mos, matrix_ks, ev_shift, eos_shift)
     765             : 
     766             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     767             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     768             :          POINTER                                         :: gs_mos
     769             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     770             :       REAL(KIND=dp), INTENT(IN)                          :: ev_shift, eos_shift
     771             : 
     772             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ev_shift_operator'
     773             : 
     774             :       INTEGER                                            :: handle, ispin, n_spins, na, nb, nhomo, &
     775             :                                                             nl, nos, nrow, nu, nvirt
     776             :       TYPE(cp_fm_struct_type), POINTER                   :: fmstruct
     777             :       TYPE(cp_fm_type)                                   :: cmos, cvec
     778             :       TYPE(cp_fm_type), POINTER                          :: coeff
     779          16 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     780             :       TYPE(dbcsr_type), POINTER                          :: smat
     781          16 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     782             : 
     783          16 :       CALL timeset(routineN, handle)
     784             : 
     785          16 :       n_spins = SIZE(gs_mos)
     786          16 :       CPASSERT(n_spins == SIZE(matrix_ks))
     787             : 
     788          16 :       IF (eos_shift /= 0.0_dp .AND. n_spins > 1) THEN
     789           0 :          CPABORT("eos_shift not implemented")
     790           0 :          CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
     791           0 :          smat => matrix_s(1)%matrix
     792           0 :          CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
     793           0 :          CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
     794           0 :          nl = MIN(na, nb)
     795           0 :          nu = MAX(na, nb)
     796             :          ! open shell orbital shift
     797           0 :          DO ispin = 1, n_spins
     798           0 :             coeff => gs_mos(ispin)%mos_occ
     799           0 :             CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
     800           0 :             IF (nhomo == nu) THEN
     801             :                ! downshift with -eos_shift using occupied orbitals
     802           0 :                nos = nu - nl
     803           0 :                CALL cp_fm_create(cmos, fmstruct)
     804           0 :                CALL cp_fm_get_info(coeff, nrow_global=nrow)
     805           0 :                CALL cp_fm_to_fm_submat(coeff, cmos, nrow, nos, 1, nl + 1, 1, 1)
     806           0 :                CALL cp_fm_create(cvec, fmstruct)
     807           0 :                CALL cp_dbcsr_sm_fm_multiply(smat, cmos, cvec, nos, 1.0_dp, 0.0_dp)
     808             :                CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
     809           0 :                                           alpha=-eos_shift, keep_sparsity=.TRUE.)
     810           0 :                CALL cp_fm_release(cmos)
     811           0 :                CALL cp_fm_release(cvec)
     812             :             ELSE
     813             :                ! upshift with eos_shift using virtual orbitals
     814           0 :                coeff => gs_mos(ispin)%mos_virt
     815           0 :                CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
     816           0 :                nos = nu - nhomo
     817           0 :                CPASSERT(nvirt >= nos)
     818           0 :                CALL cp_fm_create(cvec, fmstruct)
     819           0 :                CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
     820             :                CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
     821           0 :                                           alpha=eos_shift, keep_sparsity=.TRUE.)
     822           0 :                CALL cp_fm_release(cvec)
     823             :             END IF
     824             :          END DO
     825             :          ! virtual shift
     826           0 :          IF (ev_shift /= 0.0_dp) THEN
     827           0 :             DO ispin = 1, n_spins
     828             :                CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
     829           0 :                               alpha_scalar=1.0_dp, beta_scalar=ev_shift)
     830           0 :                coeff => gs_mos(ispin)%mos_occ
     831           0 :                CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
     832           0 :                CALL cp_fm_create(cvec, fmstruct)
     833           0 :                CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
     834             :                CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
     835           0 :                                           alpha=-ev_shift, keep_sparsity=.TRUE.)
     836           0 :                CALL cp_fm_release(cvec)
     837           0 :                IF (nhomo < nu) THEN
     838           0 :                   nos = nu - nhomo
     839           0 :                   coeff => gs_mos(ispin)%mos_virt
     840           0 :                   CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
     841           0 :                   CPASSERT(nvirt >= nos)
     842           0 :                   CALL cp_fm_create(cvec, fmstruct)
     843           0 :                   CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nos, 1.0_dp, 0.0_dp)
     844             :                   CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nos, &
     845           0 :                                              alpha=-ev_shift, keep_sparsity=.TRUE.)
     846           0 :                   CALL cp_fm_release(cvec)
     847             :                END IF
     848             :             END DO
     849             :          END IF
     850             :       ELSE
     851             :          ! virtual shift
     852          16 :          IF (ev_shift /= 0.0_dp) THEN
     853          16 :             CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
     854          16 :             smat => matrix_s(1)%matrix
     855          32 :             DO ispin = 1, n_spins
     856             :                CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
     857          16 :                               alpha_scalar=1.0_dp, beta_scalar=ev_shift)
     858          16 :                coeff => gs_mos(ispin)%mos_occ
     859          16 :                CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
     860          16 :                CALL cp_fm_create(cvec, fmstruct)
     861          16 :                CALL cp_dbcsr_sm_fm_multiply(smat, coeff, cvec, nhomo, 1.0_dp, 0.0_dp)
     862             :                CALL cp_dbcsr_plus_fm_fm_t(matrix_ks(ispin)%matrix, matrix_v=cvec, ncol=nhomo, &
     863          16 :                                           alpha=-ev_shift, keep_sparsity=.TRUE.)
     864          48 :                CALL cp_fm_release(cvec)
     865             :             END DO
     866             :          END IF
     867             :       END IF
     868             :       ! set eigenvalues
     869          16 :       IF (eos_shift == 0.0_dp .OR. n_spins == 1) THEN
     870          32 :          DO ispin = 1, n_spins
     871          32 :             IF (ALLOCATED(gs_mos(ispin)%evals_virt)) THEN
     872        1332 :                gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
     873             :             END IF
     874             :          END DO
     875             :       ELSE
     876           0 :          CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=na)
     877           0 :          CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nb)
     878           0 :          nl = MIN(na, nb)
     879           0 :          nu = MAX(na, nb)
     880           0 :          nos = nu - nl
     881           0 :          IF (na == nu) THEN
     882           0 :             IF (ALLOCATED(gs_mos(1)%evals_occ)) THEN
     883           0 :                gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
     884             :             END IF
     885           0 :             IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
     886           0 :                gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
     887             :             END IF
     888           0 :             IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
     889           0 :                gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
     890           0 :                gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
     891             :             END IF
     892             :          ELSE
     893           0 :             IF (ALLOCATED(gs_mos(1)%evals_virt)) THEN
     894           0 :                gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
     895           0 :                gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
     896             :             END IF
     897           0 :             IF (ALLOCATED(gs_mos(2)%evals_occ)) THEN
     898           0 :                gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
     899             :             END IF
     900           0 :             IF (ALLOCATED(gs_mos(2)%evals_virt)) THEN
     901           0 :                gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
     902             :             END IF
     903             :          END IF
     904             :       END IF
     905             : 
     906          16 :       CALL timestop(handle)
     907             : 
     908          16 :    END SUBROUTINE ev_shift_operator
     909             : 
     910             : ! **************************************************************************************************
     911             : !> \brief Generate missed guess vectors.
     912             : !> \param evects   guess vectors distributed across all processors (initialised on exit)
     913             : !> \param evals    guessed transition energies (initialised on exit)
     914             : !> \param gs_mos   occupied and virtual molecular orbitals optimised for the ground state
     915             : !> \param log_unit output unit
     916             : !> \par History
     917             : !>    * 05.2016 created as tddfpt_guess() [Sergey Chulkov]
     918             : !>    * 06.2016 renamed, altered prototype, supports spin-polarised density [Sergey Chulkov]
     919             : !>    * 01.2017 simplified prototype, do not compute all possible singly-excited states
     920             : !>              [Sergey Chulkov]
     921             : !> \note \parblock
     922             : !>       Based on the subroutine co_initial_guess() which was originally created by
     923             : !>       Thomas Chassaing on 06.2003.
     924             : !>
     925             : !>       Only not associated guess vectors 'evects(spin, state)%matrix' are allocated and
     926             : !>       initialised; associated vectors assumed to be initialised elsewhere (e.g. using
     927             : !>       a restart file).
     928             : !>       \endparblock
     929             : ! **************************************************************************************************
     930        1062 :    SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit)
     931             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout)   :: evects
     932             :       REAL(kind=dp), DIMENSION(:), INTENT(inout)         :: evals
     933             :       TYPE(tddfpt_ground_state_mos), &
     934             :          DIMENSION(SIZE(evects, 1)), INTENT(in)          :: gs_mos
     935             :       INTEGER, INTENT(in)                                :: log_unit
     936             : 
     937             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_guess_vectors'
     938             : 
     939             :       CHARACTER(len=5)                                   :: spin_label
     940             :       INTEGER                                            :: handle, imo_occ, imo_virt, ind, ispin, &
     941             :                                                             istate, jspin, nspins, nstates, &
     942             :                                                             nstates_occ_virt_alpha, &
     943             :                                                             nstates_selected
     944        1062 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: inds
     945             :       INTEGER, DIMENSION(maxspins)                       :: nmo_occ_avail, nmo_occ_selected, &
     946             :                                                             nmo_virt_selected
     947             :       REAL(kind=dp)                                      :: e_occ
     948        1062 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: e_virt_minus_occ
     949        3186 :       TYPE(cp_fm_struct_p_type), DIMENSION(maxspins)     :: fm_struct_evects
     950             : 
     951        1062 :       CALL timeset(routineN, handle)
     952             : 
     953        1062 :       nspins = SIZE(evects, 1)
     954        1062 :       nstates = SIZE(evects, 2)
     955             : 
     956             :       IF (debug_this_module) THEN
     957             :          CPASSERT(nstates > 0)
     958             :          CPASSERT(nspins == 1 .OR. nspins == 2)
     959             :       END IF
     960             : 
     961        2248 :       DO ispin = 1, nspins
     962             :          ! number of occupied orbitals for each spin component
     963        1186 :          nmo_occ_avail(ispin) = SIZE(gs_mos(ispin)%evals_occ)
     964             :          ! number of occupied and virtual orbitals which can potentially
     965             :          ! contribute to the excited states in question.
     966        1186 :          nmo_occ_selected(ispin) = MIN(nmo_occ_avail(ispin), nstates)
     967        1186 :          nmo_virt_selected(ispin) = MIN(SIZE(gs_mos(ispin)%evals_virt), nstates)
     968             : 
     969        2248 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct_evects(ispin)%struct)
     970             :       END DO
     971             : 
     972             :       ! TO DO: the variable 'nstates_selected' should probably be declared as INTEGER(kind=int_8),
     973             :       !        however we need a special version of the subroutine sort() in order to do so
     974        2248 :       nstates_selected = DOT_PRODUCT(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
     975             : 
     976        3186 :       ALLOCATE (inds(nstates_selected))
     977        3186 :       ALLOCATE (e_virt_minus_occ(nstates_selected))
     978             : 
     979        1062 :       istate = 0
     980        2248 :       DO ispin = 1, nspins
     981        5066 :          DO imo_occ = 1, nmo_occ_selected(ispin)
     982             :             ! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element)
     983        2818 :             e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1)
     984             : 
     985       13930 :             DO imo_virt = 1, nmo_virt_selected(ispin)
     986        9926 :                istate = istate + 1
     987       12744 :                e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ
     988             :             END DO
     989             :          END DO
     990             :       END DO
     991             : 
     992             :       IF (debug_this_module) THEN
     993             :          CPASSERT(istate == nstates_selected)
     994             :       END IF
     995             : 
     996        1062 :       CALL sort(e_virt_minus_occ, nstates_selected, inds)
     997             : 
     998        1062 :       IF (nspins == 1) THEN
     999         938 :          ispin = 1
    1000         938 :          spin_label = '     '
    1001             :       END IF
    1002             : 
    1003        1062 :       nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
    1004        1062 :       IF (log_unit > 0) THEN
    1005         531 :          WRITE (log_unit, "(1X,A)") "", &
    1006         531 :             "-------------------------------------------------------------------------------", &
    1007         531 :             "-                            TDDFPT Initial Guess                             -", &
    1008        1062 :             "-------------------------------------------------------------------------------"
    1009         531 :          WRITE (log_unit, '(T11,A)') "State         Occupied      ->      Virtual          Excitation"
    1010         531 :          WRITE (log_unit, '(T11,A)') "number         orbital              orbital          energy (eV)"
    1011         531 :          WRITE (log_unit, '(1X,79("-"))')
    1012             :       END IF
    1013             : 
    1014        3936 :       DO istate = 1, nstates
    1015        3936 :          IF (ASSOCIATED(evects(1, istate)%matrix_struct)) THEN
    1016           6 :             IF (log_unit > 0) &
    1017             :                WRITE (log_unit, '(T7,I8,T28,A19,T60,F14.5)') &
    1018           3 :                istate, "***  restarted  ***", evals(istate)*evolt
    1019             :          ELSE
    1020        2868 :             ind = inds(istate) - 1
    1021        2868 :             IF (nspins > 1) THEN
    1022         372 :                IF (ind < nstates_occ_virt_alpha) THEN
    1023         152 :                   ispin = 1
    1024         152 :                   spin_label = '(alp)'
    1025             :                ELSE
    1026         220 :                   ispin = 2
    1027         220 :                   ind = ind - nstates_occ_virt_alpha
    1028         220 :                   spin_label = '(bet)'
    1029             :                END IF
    1030             :             END IF
    1031             : 
    1032        2868 :             imo_occ = nmo_occ_avail(ispin) - ind/nmo_virt_selected(ispin)
    1033        2868 :             imo_virt = MOD(ind, nmo_virt_selected(ispin)) + 1
    1034        2868 :             evals(istate) = e_virt_minus_occ(istate)
    1035             : 
    1036        2868 :             IF (log_unit > 0) &
    1037             :                WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
    1038        1434 :                istate, imo_occ, spin_label, nmo_occ_avail(ispin) + imo_virt, spin_label, e_virt_minus_occ(istate)*evolt
    1039             : 
    1040        6108 :             DO jspin = 1, nspins
    1041             :                ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
    1042        3240 :                CALL cp_fm_create(evects(jspin, istate), fm_struct_evects(jspin)%struct)
    1043        3240 :                CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
    1044             : 
    1045        3240 :                IF (jspin == ispin) &
    1046             :                   CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, evects(ispin, istate), &
    1047        5736 :                                    ncol=1, source_start=imo_virt, target_start=imo_occ)
    1048             :             END DO
    1049             :          END IF
    1050             :       END DO
    1051             : 
    1052        1062 :       IF (log_unit > 0) THEN
    1053         531 :          WRITE (log_unit, '(/,T7,A,T50,I24)') 'Number of active states:', tddfpt_total_number_of_states(gs_mos)
    1054             :          WRITE (log_unit, "(1X,A)") &
    1055         531 :             "-------------------------------------------------------------------------------"
    1056             :       END IF
    1057             : 
    1058        1062 :       DEALLOCATE (e_virt_minus_occ)
    1059        1062 :       DEALLOCATE (inds)
    1060             : 
    1061        1062 :       CALL timestop(handle)
    1062             : 
    1063        1062 :    END SUBROUTINE tddfpt_guess_vectors
    1064             : 
    1065             : END MODULE qs_tddfpt2_utils

Generated by: LCOV version 1.15