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

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : MODULE qs_tddfpt2_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_api,                    ONLY: dbcsr_add,&
      15              :                                               dbcsr_copy,&
      16              :                                               dbcsr_get_info,&
      17              :                                               dbcsr_init_p,&
      18              :                                               dbcsr_p_type,&
      19              :                                               dbcsr_type
      20              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      21              :                                               cp_dbcsr_plus_fm_fm_t,&
      22              :                                               cp_dbcsr_sm_fm_multiply,&
      23              :                                               dbcsr_allocate_matrix_set
      24              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_triangular_invert
      25              :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose
      26              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      27              :                                               cp_fm_struct_p_type,&
      28              :                                               cp_fm_struct_release,&
      29              :                                               cp_fm_struct_type
      30              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      31              :                                               cp_fm_get_info,&
      32              :                                               cp_fm_release,&
      33              :                                               cp_fm_set_all,&
      34              :                                               cp_fm_to_fm,&
      35              :                                               cp_fm_to_fm_submat,&
      36              :                                               cp_fm_type
      37              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      38              :                                               cp_logger_get_default_io_unit,&
      39              :                                               cp_logger_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_ks_matrix,                   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         1076 :    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         1076 :       REAL(kind=dp), DIMENSION(:), POINTER               :: evals_virt_spin
     112              :       TYPE(cell_type), POINTER                           :: cell
     113         1076 :       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         1076 :          TARGET                                          :: mos_virt
     117              :       TYPE(cp_fm_type), POINTER                          :: mos_virt_spin
     118         1076 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s
     119              :       TYPE(dft_control_type), POINTER                    :: dft_control
     120         1076 :       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         1076 :       CALL timeset(routineN, handle)
     126              : 
     127              :       CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
     128         1076 :                       matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
     129         1076 :       tddfpt_control => dft_control%tddfpt2_control
     130              : 
     131         1076 :       CPASSERT(.NOT. ASSOCIATED(gs_mos))
     132              :       ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
     133         1076 :       nspins = dft_control%nspins
     134         4428 :       ALLOCATE (gs_mos(nspins))
     135              : 
     136              :       ! check if virtuals should be constructed for NAMD interface with NEWTONX
     137         1076 :       print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
     138         1076 :       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         1076 :       NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
     144         1076 :       IF (ASSOCIATED(scf_env)) THEN
     145         1076 :          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         2276 :       DO ispin = 1, nspins
     168         1200 :          IF (ASSOCIATED(evals_virt)) THEN
     169            0 :             evals_virt_spin => evals_virt(ispin)%array
     170              :          ELSE
     171         1200 :             NULLIFY (evals_virt_spin)
     172              :          END IF
     173         1200 :          IF (ALLOCATED(mos_virt)) THEN
     174            0 :             mos_virt_spin => mos_virt(ispin)
     175              :          ELSE
     176         1200 :             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         2276 :                                            qs_env=qs_env)
     184              :       END DO
     185              : 
     186         1076 :       moc = 0
     187         1076 :       mvt = 0
     188         2276 :       DO ispin = 1, nspins
     189         1200 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
     190         2276 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
     191              :       END DO
     192         1076 :       IF (iounit > 0) THEN
     193          538 :          WRITE (iounit, "(T2,A,T36,A)") "TDDFPT| Molecular Orbitals:", &
     194         1076 :             " Spin       AOs       Occ      Virt     Total"
     195         1138 :          DO ispin = 1, nspins
     196          600 :             WRITE (iounit, "(T2,A,T37,I4,4I10)") "TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
     197         1738 :                mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
     198              :          END DO
     199              :       END IF
     200              : 
     201         1076 :       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         1076 :       CALL cp_fm_release(mos_virt)
     209              : 
     210         1076 :       CALL timestop(handle)
     211              : 
     212         3228 :    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         1200 :    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         1200 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: minrow_neg_array, minrow_pos_array, &
     253         1200 :                                                             sum_sign_array
     254         1200 :       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         1200 :          POINTER                                         :: my_block
     259         1200 :       REAL(kind=dp), DIMENSION(:), POINTER               :: mo_evals_extended, mo_occ_extended, &
     260         1200 :                                                             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         1200 :       CALL timeset(routineN, handle)
     272              : 
     273         1200 :       NULLIFY (logger)
     274         1200 :       logger => cp_get_default_logger()
     275         1200 :       iounit = cp_logger_get_default_io_unit(logger)
     276              : 
     277         1200 :       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         1200 :                       nelectron=nelectrons, occupation_numbers=mo_occ_scf)
     281              : 
     282         1200 :       print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
     283         1200 :       CALL section_vals_val_get(print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
     284              : 
     285         1200 :       nmo_virt = nao - nmo_occ
     286         1200 :       IF (nlumo >= 0) &
     287            0 :          nmo_virt = MIN(nmo_virt, nlumo)
     288              : 
     289         1200 :       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         1200 :       do_eigen = .FALSE.
     294              :       ! diagonalise the Kohn-Sham matrix one more time if the number of available unoccupied states are too small
     295         1200 :       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         1200 :          IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .TRUE.
     300              :       END IF
     301              : 
     302              :       ! ++ allocate storage space for gs_mos
     303         1200 :       NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
     304         1200 :       CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
     305         1200 :       CALL cp_fm_struct_create(ao_mo_virt_fm_struct, nrow_global=nao, ncol_global=nmo_virt, context=blacs_env)
     306              : 
     307         1200 :       NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
     308         1200 :       ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
     309         1200 :       CALL cp_fm_create(gs_mos%mos_occ, ao_mo_occ_fm_struct)
     310         1200 :       CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
     311              : 
     312         3600 :       ALLOCATE (gs_mos%evals_occ(nmo_occ))
     313         3600 :       ALLOCATE (gs_mos%evals_virt(nmo_virt))
     314         2400 :       ALLOCATE (gs_mos%phases_occ(nmo_occ))
     315         2400 :       ALLOCATE (gs_mos%phases_virt(nmo_virt))
     316              : 
     317              :       ! ++ nullify pointers
     318         1200 :       NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
     319         1200 :       NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
     320         1200 :       CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
     321              : 
     322         1200 :       IF (do_eigen) THEN
     323              :          ! ++ set of molecular orbitals
     324         1200 :          CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
     325         1200 :          ALLOCATE (mos_extended)
     326              :          CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
     327         1200 :                               REAL(nelectrons, dp), maxocc, flexible_electron_count=0.0_dp)
     328         1200 :          CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name="mos-extended")
     329         1200 :          CALL cp_fm_struct_release(wfn_fm_struct)
     330              :          CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
     331         1200 :                          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         7596 :          DO imo = 1, nmo_scf
     339         7596 :             mo_occ_extended(imo) = mo_occ_scf(imo)
     340              :          END DO
     341        20850 :          mo_occ_extended(nmo_scf + 1:) = 0.0_dp
     342              : 
     343              :          ! ++ allocate temporary matrices
     344         1200 :          CALL cp_fm_create(matrix_ks_fm, ao_ao_fm_struct)
     345         1200 :          CALL cp_fm_create(ortho_fm, ao_ao_fm_struct)
     346         1200 :          CALL cp_fm_create(work_fm, ao_ao_fm_struct)
     347         1200 :          CALL cp_fm_struct_release(ao_ao_fm_struct)
     348              : 
     349              :          ! some stuff from the subroutine general_eigenproblem()
     350         1200 :          CALL copy_dbcsr_to_fm(matrix_s, ortho_fm)
     351         1200 :          CALL copy_dbcsr_to_fm(matrix_ks, matrix_ks_fm)
     352              : 
     353         1200 :          IF (cholesky_method == cholesky_dbcsr) THEN
     354            0 :             CPABORT('CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
     355         1200 :          ELSE IF (cholesky_method == cholesky_off) THEN
     356            0 :             CPABORT('CHOLESKY OFF is not implemented in TDDFT.')
     357              :          ELSE
     358         1200 :             CALL cp_fm_cholesky_decompose(ortho_fm)
     359         1200 :             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         1200 :             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         1200 :                              do_level_shift=.FALSE., level_shift=0.0_dp, use_jacobi=.FALSE.)
     369              :          END IF
     370              : 
     371              :          ! -- clean up needless matrices
     372         1200 :          CALL cp_fm_release(work_fm)
     373         1200 :          CALL cp_fm_release(ortho_fm)
     374         1200 :          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         1200 :       CALL cp_fm_create(work_fm, ao_mo_occ_fm_struct)
     384         1200 :       CALL cp_fm_struct_release(ao_mo_occ_fm_struct)
     385              : 
     386         1200 :       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         1200 :                           row_indices=row_indices, col_indices=col_indices, local_data=my_block)
     389              : 
     390         6000 :       ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
     391         7576 :       minrow_neg_array(:) = nao
     392         7576 :       minrow_pos_array(:) = nao
     393         7576 :       sum_sign_array(:) = 0
     394         7576 :       DO icol_local = 1, ncol_local
     395         6376 :          icol_global = col_indices(icol_local)
     396              : 
     397       224116 :          DO irow_local = 1, nrow_local
     398       216540 :             element = my_block(irow_local, icol_local)
     399              : 
     400       216540 :             sign_int = 0
     401       216540 :             IF (element >= eps_dp) THEN
     402              :                sign_int = 1
     403       109407 :             ELSE IF (element <= -eps_dp) THEN
     404       108877 :                sign_int = -1
     405              :             END IF
     406              : 
     407       216540 :             sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
     408              : 
     409       216540 :             irow_global = row_indices(irow_local)
     410       222916 :             IF (sign_int > 0) THEN
     411       107133 :                IF (minrow_pos_array(icol_global) > irow_global) &
     412         6259 :                   minrow_pos_array(icol_global) = irow_global
     413       109407 :             ELSE IF (sign_int < 0) THEN
     414       108877 :                IF (minrow_neg_array(icol_global) > irow_global) &
     415         6132 :                   minrow_neg_array(icol_global) = irow_global
     416              :             END IF
     417              :          END DO
     418              :       END DO
     419              : 
     420         1200 :       CALL para_env%sum(sum_sign_array)
     421         1200 :       CALL para_env%min(minrow_neg_array)
     422         1200 :       CALL para_env%min(minrow_pos_array)
     423              : 
     424         7576 :       DO icol_local = 1, nmo_occ
     425         7576 :          IF (sum_sign_array(icol_local) > 0) THEN
     426              :             ! most of the expansion coefficients are positive => MO's phase = +1
     427         2830 :             gs_mos%phases_occ(icol_local) = 1.0_dp
     428         3546 :          ELSE IF (sum_sign_array(icol_local) < 0) THEN
     429              :             ! most of the expansion coefficients are negative => MO's phase = -1
     430         3302 :             gs_mos%phases_occ(icol_local) = -1.0_dp
     431              :          ELSE
     432              :             ! equal number of positive and negative expansion coefficients
     433          244 :             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          134 :                gs_mos%phases_occ(icol_local) = 1.0_dp
     437              :             ELSE
     438              :                ! MO's phase = -1
     439          110 :                gs_mos%phases_occ(icol_local) = -1.0_dp
     440              :             END IF
     441              :          END IF
     442              :       END DO
     443              : 
     444         1200 :       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         1200 :       CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
     448         7576 :       gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
     449              : 
     450         1200 :       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         1200 :          CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
     459        20870 :          gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
     460              :       END IF
     461              : 
     462         1200 :       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         1200 :       CALL cp_fm_struct_release(ao_mo_virt_fm_struct) ! here after print_phases
     528              : 
     529         1200 :       CALL cp_fm_release(work_fm)
     530              : 
     531         1200 :       IF (do_eigen) THEN
     532         1200 :          CALL deallocate_mo_set(mos_extended)
     533         1200 :          DEALLOCATE (mos_extended)
     534              :       END IF
     535              : 
     536         1200 :       CALL timestop(handle)
     537              : 
     538         7200 :    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         1200 :    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         1200 :       CALL timeset(routineN, handle)
     554              : 
     555         1200 :       IF (ALLOCATED(gs_mos%phases_occ)) &
     556         1200 :          DEALLOCATE (gs_mos%phases_occ)
     557              : 
     558         1200 :       IF (ALLOCATED(gs_mos%evals_virt)) &
     559         1200 :          DEALLOCATE (gs_mos%evals_virt)
     560              : 
     561         1200 :       IF (ALLOCATED(gs_mos%evals_occ)) &
     562         1200 :          DEALLOCATE (gs_mos%evals_occ)
     563              : 
     564         1200 :       IF (ALLOCATED(gs_mos%phases_virt)) &
     565         1200 :          DEALLOCATE (gs_mos%phases_virt)
     566              : 
     567         1200 :       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         1200 :       IF (ASSOCIATED(gs_mos%mos_virt)) THEN
     573         1200 :          CALL cp_fm_release(gs_mos%mos_virt)
     574         1200 :          DEALLOCATE (gs_mos%mos_virt)
     575              :       END IF
     576              : 
     577         1200 :       IF (ASSOCIATED(gs_mos%mos_occ)) THEN
     578         1200 :          CALL cp_fm_release(gs_mos%mos_occ)
     579         1200 :          DEALLOCATE (gs_mos%mos_occ)
     580              :       END IF
     581              : 
     582         1200 :       CALL timestop(handle)
     583         1200 :    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         1076 :    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         1076 :       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         1076 :       CALL timeset(routineN, handle)
     614              : 
     615         1076 :       NULLIFY (logger)
     616         1076 :       logger => cp_get_default_logger()
     617         1076 :       iounit = cp_logger_get_default_io_unit(logger)
     618              : 
     619         1076 :       CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
     620         1076 :       tddfpt_control => dft_control%tddfpt2_control
     621              : 
     622              :       ! obtain corrected KS-matrix
     623              :       ! We should 'save' the energy values?
     624         1076 :       nspins = SIZE(gs_mos)
     625         1076 :       NULLIFY (matrix_ks_oep)
     626         1076 :       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          102 :             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         1076 :       CALL timestop(handle)
     724              : 
     725         1076 :    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         1708 :    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         1708 :       nstates_total = 0
     742         1708 :       nspins = SIZE(gs_mos)
     743              : 
     744         3602 :       DO ispin = 1, nspins
     745              :          nstates_total = nstates_total + &
     746              :                          SIZE(gs_mos(ispin)%evals_occ, kind=int_8)* &
     747         3602 :                          SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
     748              :       END DO
     749         1708 :    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         1084 :    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         1084 :       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         1084 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: e_virt_minus_occ
     949         3252 :       TYPE(cp_fm_struct_p_type), DIMENSION(maxspins)     :: fm_struct_evects
     950              : 
     951         1084 :       CALL timeset(routineN, handle)
     952              : 
     953         1084 :       nspins = SIZE(evects, 1)
     954         1084 :       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         2292 :       DO ispin = 1, nspins
     962              :          ! number of occupied orbitals for each spin component
     963         1208 :          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         1208 :          nmo_occ_selected(ispin) = MIN(nmo_occ_avail(ispin), nstates)
     967         1208 :          nmo_virt_selected(ispin) = MIN(SIZE(gs_mos(ispin)%evals_virt), nstates)
     968              : 
     969         2292 :          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         2292 :       nstates_selected = DOT_PRODUCT(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
     975              : 
     976         3252 :       ALLOCATE (inds(nstates_selected))
     977         3252 :       ALLOCATE (e_virt_minus_occ(nstates_selected))
     978              : 
     979         1084 :       istate = 0
     980         2292 :       DO ispin = 1, nspins
     981         5164 :          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         2872 :             e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1)
     984              : 
     985        14172 :             DO imo_virt = 1, nmo_virt_selected(ispin)
     986        10092 :                istate = istate + 1
     987        12964 :                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         1084 :       CALL sort(e_virt_minus_occ, nstates_selected, inds)
     997              : 
     998         1084 :       IF (nspins == 1) THEN
     999          960 :          ispin = 1
    1000          960 :          spin_label = '     '
    1001              :       END IF
    1002              : 
    1003         1084 :       nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
    1004         1084 :       IF (log_unit > 0) THEN
    1005          542 :          WRITE (log_unit, "(1X,A)") "", &
    1006          542 :             "-------------------------------------------------------------------------------", &
    1007          542 :             "-                            TDDFPT Initial Guess                             -", &
    1008         1084 :             "-------------------------------------------------------------------------------"
    1009          542 :          WRITE (log_unit, '(T11,A)') "State         Occupied      ->      Virtual          Excitation"
    1010          542 :          WRITE (log_unit, '(T11,A)') "number         orbital              orbital          energy (eV)"
    1011          542 :          WRITE (log_unit, '(1X,79("-"))')
    1012              :       END IF
    1013              : 
    1014         4018 :       DO istate = 1, nstates
    1015         4018 :          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         2928 :             ind = inds(istate) - 1
    1021         2928 :             IF (nspins > 1) THEN
    1022          372 :                IF (ind < nstates_occ_virt_alpha) THEN
    1023          158 :                   ispin = 1
    1024          158 :                   spin_label = '(alp)'
    1025              :                ELSE
    1026          214 :                   ispin = 2
    1027          214 :                   ind = ind - nstates_occ_virt_alpha
    1028          214 :                   spin_label = '(bet)'
    1029              :                END IF
    1030              :             END IF
    1031              : 
    1032         2928 :             imo_occ = nmo_occ_avail(ispin) - ind/nmo_virt_selected(ispin)
    1033         2928 :             imo_virt = MOD(ind, nmo_virt_selected(ispin)) + 1
    1034         2928 :             evals(istate) = e_virt_minus_occ(istate)
    1035              : 
    1036         2928 :             IF (log_unit > 0) &
    1037              :                WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
    1038         1464 :                istate, imo_occ, spin_label, nmo_occ_avail(ispin) + imo_virt, spin_label, e_virt_minus_occ(istate)*evolt
    1039              : 
    1040         6228 :             DO jspin = 1, nspins
    1041              :                ! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
    1042         3300 :                CALL cp_fm_create(evects(jspin, istate), fm_struct_evects(jspin)%struct)
    1043         3300 :                CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
    1044              : 
    1045         3300 :                IF (jspin == ispin) &
    1046              :                   CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, evects(ispin, istate), &
    1047         5856 :                                    ncol=1, source_start=imo_virt, target_start=imo_occ)
    1048              :             END DO
    1049              :          END IF
    1050              :       END DO
    1051              : 
    1052         1084 :       IF (log_unit > 0) THEN
    1053          542 :          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          542 :             "-------------------------------------------------------------------------------"
    1056              :       END IF
    1057              : 
    1058         1084 :       DEALLOCATE (e_virt_minus_occ)
    1059         1084 :       DEALLOCATE (inds)
    1060              : 
    1061         1084 :       CALL timestop(handle)
    1062              : 
    1063         1084 :    END SUBROUTINE tddfpt_guess_vectors
    1064              : 
    1065              : END MODULE qs_tddfpt2_utils
        

Generated by: LCOV version 2.0-1