LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:4dc10b3) Lines: 484 528 91.7 %
Date: 2024-11-21 06:45:46 Functions: 5 5 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_methods
       9             :    USE admm_methods,                    ONLY: admm_fit_mo_coeffs
      10             :    USE admm_types,                      ONLY: admm_type,&
      11             :                                               get_admm_env
      12             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      13             :    USE bibliography,                    ONLY: Grimme2013,&
      14             :                                               Grimme2016,&
      15             :                                               Iannuzzi2005,&
      16             :                                               cite_reference
      17             :    USE cell_types,                      ONLY: cell_type
      18             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      19             :    USE cp_control_types,                ONLY: dft_control_type,&
      20             :                                               tddfpt2_control_type
      21             :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      22             :    USE cp_dbcsr_operations,             ONLY: dbcsr_deallocate_matrix_set
      23             :    USE cp_fm_pool_types,                ONLY: fm_pool_create_fm
      24             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_type
      25             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      26             :                                               cp_fm_get_info,&
      27             :                                               cp_fm_release,&
      28             :                                               cp_fm_to_fm,&
      29             :                                               cp_fm_type
      30             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      31             :                                               cp_logger_get_default_io_unit,&
      32             :                                               cp_logger_type
      33             :    USE cp_output_handling,              ONLY: cp_add_iter_level,&
      34             :                                               cp_iterate,&
      35             :                                               cp_print_key_finished_output,&
      36             :                                               cp_print_key_unit_nr,&
      37             :                                               cp_rm_iter_level
      38             :    USE exstates_types,                  ONLY: excited_energy_type
      39             :    USE header,                          ONLY: tddfpt_header,&
      40             :                                               tddfpt_soc_header
      41             :    USE hfx_admm_utils,                  ONLY: aux_admm_init
      42             :    USE hfx_types,                       ONLY: compare_hfx_sections,&
      43             :                                               hfx_create
      44             :    USE input_constants,                 ONLY: &
      45             :         do_admm_aux_exch_func_none, do_admm_basis_projection, do_admm_exch_scaling_none, &
      46             :         do_admm_purify_none, do_potential_truncated, tddfpt_dipole_velocity, tddfpt_kernel_full, &
      47             :         tddfpt_kernel_none, tddfpt_kernel_stda
      48             :    USE input_section_types,             ONLY: section_vals_get,&
      49             :                                               section_vals_get_subs_vals,&
      50             :                                               section_vals_type,&
      51             :                                               section_vals_val_get,&
      52             :                                               section_vals_val_set
      53             :    USE kinds,                           ONLY: dp
      54             :    USE lri_environment_methods,         ONLY: lri_print_stat
      55             :    USE lri_environment_types,           ONLY: lri_density_release,&
      56             :                                               lri_env_release
      57             :    USE machine,                         ONLY: m_flush
      58             :    USE message_passing,                 ONLY: mp_para_env_type
      59             :    USE min_basis_set,                   ONLY: create_minbas_set
      60             :    USE particle_types,                  ONLY: particle_type
      61             :    USE physcon,                         ONLY: evolt
      62             :    USE qs_environment_types,            ONLY: get_qs_env,&
      63             :                                               qs_environment_type
      64             :    USE qs_kernel_methods,               ONLY: create_fxc_kernel,&
      65             :                                               create_kernel_env
      66             :    USE qs_kernel_types,                 ONLY: full_kernel_env_type,&
      67             :                                               kernel_env_type,&
      68             :                                               release_kernel_env
      69             :    USE qs_kind_types,                   ONLY: qs_kind_type
      70             :    USE qs_mo_types,                     ONLY: mo_set_type
      71             :    USE qs_scf_methods,                  ONLY: eigensolver
      72             :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      73             :    USE qs_tddfpt2_assign,               ONLY: assign_state
      74             :    USE qs_tddfpt2_densities,            ONLY: tddfpt_construct_aux_fit_density,&
      75             :                                               tddfpt_construct_ground_state_orb_density
      76             :    USE qs_tddfpt2_eigensolver,          ONLY: tddfpt_davidson_solver,&
      77             :                                               tddfpt_orthogonalize_psi1_psi0,&
      78             :                                               tddfpt_orthonormalize_psi1_psi1
      79             :    USE qs_tddfpt2_forces,               ONLY: tddfpt_forces_main
      80             :    USE qs_tddfpt2_fprint,               ONLY: tddfpt_print_forces
      81             :    USE qs_tddfpt2_lri_utils,            ONLY: tddfpt2_lri_init
      82             :    USE qs_tddfpt2_properties,           ONLY: tddfpt_dipole_operator,&
      83             :                                               tddfpt_print_excitation_analysis,&
      84             :                                               tddfpt_print_nto_analysis,&
      85             :                                               tddfpt_print_summary
      86             :    USE qs_tddfpt2_restart,              ONLY: tddfpt_read_restart,&
      87             :                                               tddfpt_write_newtonx_output,&
      88             :                                               tddfpt_write_restart
      89             :    USE qs_tddfpt2_soc,                  ONLY: tddfpt_soc
      90             :    USE qs_tddfpt2_stda_types,           ONLY: allocate_stda_env,&
      91             :                                               deallocate_stda_env,&
      92             :                                               stda_env_type,&
      93             :                                               stda_init_param
      94             :    USE qs_tddfpt2_stda_utils,           ONLY: get_lowdin_mo_coefficients,&
      95             :                                               stda_init_matrices
      96             :    USE qs_tddfpt2_subgroups,            ONLY: tddfpt_sub_env_init,&
      97             :                                               tddfpt_sub_env_release,&
      98             :                                               tddfpt_subgroup_env_type
      99             :    USE qs_tddfpt2_types,                ONLY: hfxsr_create_work_matrices,&
     100             :                                               stda_create_work_matrices,&
     101             :                                               tddfpt_create_work_matrices,&
     102             :                                               tddfpt_ground_state_mos,&
     103             :                                               tddfpt_release_work_matrices,&
     104             :                                               tddfpt_work_matrices
     105             :    USE qs_tddfpt2_utils,                ONLY: tddfpt_guess_vectors,&
     106             :                                               tddfpt_init_mos,&
     107             :                                               tddfpt_oecorr,&
     108             :                                               tddfpt_release_ground_state_mos
     109             :    USE string_utilities,                ONLY: integer_to_string
     110             :    USE xc_write_output,                 ONLY: xc_write
     111             : #include "./base/base_uses.f90"
     112             : 
     113             :    IMPLICIT NONE
     114             : 
     115             :    PRIVATE
     116             : 
     117             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_methods'
     118             : 
     119             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
     120             :    ! number of first derivative components (3: d/dx, d/dy, d/dz)
     121             :    INTEGER, PARAMETER, PRIVATE          :: nderivs = 3
     122             :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
     123             : 
     124             :    PUBLIC :: tddfpt, tddfpt_energies, tddfpt_input
     125             : 
     126             : ! **************************************************************************************************
     127             : 
     128             : CONTAINS
     129             : 
     130             : ! **************************************************************************************************
     131             : !> \brief Perform TDDFPT calculation.
     132             : !> \param qs_env  Quickstep environment
     133             : !> \param calc_forces ...
     134             : !> \par History
     135             : !>    * 05.2016 created [Sergey Chulkov]
     136             : !>    * 06.2016 refactored to be used with Davidson eigensolver [Sergey Chulkov]
     137             : !>    * 03.2017 cleaned and refactored [Sergey Chulkov]
     138             : !> \note Based on the subroutines tddfpt_env_init(), and tddfpt_env_deallocate().
     139             : ! **************************************************************************************************
     140        1054 :    SUBROUTINE tddfpt(qs_env, calc_forces)
     141             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     142             :       LOGICAL, INTENT(IN)                                :: calc_forces
     143             : 
     144             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt'
     145             : 
     146             :       INTEGER                                            :: handle, ispin, istate, log_unit, mult, &
     147             :                                                             my_state, nao, nocc, nspins, &
     148             :                                                             nstate_max, nstates, nvirt, old_state
     149             :       INTEGER, DIMENSION(maxspins)                       :: nactive
     150             :       LOGICAL                                            :: do_admm, do_exck, do_hfx, do_hfxlr, &
     151             :                                                             do_hfxsr, do_soc, lmult_tmp, &
     152             :                                                             state_change
     153             :       REAL(kind=dp)                                      :: gsmin, gsval, xsval
     154        1054 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: evals, ostrength
     155        1054 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     156             :       TYPE(cell_type), POINTER                           :: cell
     157             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     158             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     159        1054 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: my_mos
     160        1054 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: dipole_op_mos_occ, evects, S_evects
     161             :       TYPE(cp_logger_type), POINTER                      :: logger
     162        1054 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_ks_oep, matrix_s, &
     163        1054 :                                                             matrix_s_aux_fit, &
     164        1054 :                                                             matrix_s_aux_fit_vs_orb
     165             :       TYPE(dft_control_type), POINTER                    :: dft_control
     166             :       TYPE(excited_energy_type), POINTER                 :: ex_env
     167             :       TYPE(full_kernel_env_type), TARGET                 :: full_kernel_env, kernel_env_admm_aux
     168             :       TYPE(kernel_env_type)                              :: kernel_env
     169        1054 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos, mos_aux_fit
     170             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     171        1054 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     172        1054 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     173             :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     174             :       TYPE(section_vals_type), POINTER                   :: hfxsr_section, kernel_section, &
     175             :                                                             lri_section, soc_section, &
     176             :                                                             tddfpt_print_section, tddfpt_section, &
     177             :                                                             xc_section
     178             :       TYPE(stda_env_type), TARGET                        :: stda_kernel
     179             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     180             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     181        1054 :          POINTER                                         :: gs_mos
     182        1054 :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
     183        1054 :       TYPE(tddfpt_work_matrices)                         :: work_matrices
     184             : 
     185        1054 :       CALL timeset(routineN, handle)
     186             : 
     187        1054 :       NULLIFY (logger)
     188        1054 :       logger => cp_get_default_logger()
     189             : 
     190             :       ! input section print/xc
     191        1054 :       NULLIFY (tddfpt_section)
     192        1054 :       tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")
     193             : 
     194             :       CALL tddfpt_input(qs_env, do_hfx, do_admm, do_exck, &
     195             :                         do_hfxsr, do_hfxlr, xc_section, tddfpt_print_section, &
     196        1054 :                         lri_section, hfxsr_section)
     197             : 
     198             :       log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, "PROGRAM_BANNER", &
     199        1054 :                                       extension=".tddfptLog")
     200             : 
     201             :       CALL get_qs_env(qs_env, &
     202             :                       blacs_env=blacs_env, &
     203             :                       cell=cell, &
     204             :                       dft_control=dft_control, &
     205             :                       matrix_ks=matrix_ks, &
     206             :                       matrix_s=matrix_s, &
     207             :                       mos=mos, &
     208        1054 :                       scf_env=scf_env)
     209        1054 :       tddfpt_control => dft_control%tddfpt2_control
     210        1054 :       tddfpt_control%do_hfx = do_hfx
     211        1054 :       tddfpt_control%do_admm = do_admm
     212        1054 :       tddfpt_control%do_hfxsr = do_hfxsr
     213        1054 :       tddfpt_control%hfxsr_primbas = 0
     214        1054 :       tddfpt_control%hfxsr_re_int = .TRUE.
     215        1054 :       tddfpt_control%do_hfxlr = do_hfxlr
     216        1054 :       tddfpt_control%do_exck = do_exck
     217        1054 :       IF (tddfpt_control%do_hfxlr) THEN
     218           6 :          kernel_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL%HFXLR")
     219           6 :          CALL section_vals_val_get(kernel_section, "RCUT", r_val=tddfpt_control%hfxlr_rcut)
     220           6 :          CALL section_vals_val_get(kernel_section, "SCALE", r_val=tddfpt_control%hfxlr_scale)
     221             :       END IF
     222             : 
     223        1054 :       soc_section => section_vals_get_subs_vals(tddfpt_section, "SOC")
     224        1054 :       CALL section_vals_get(soc_section, explicit=do_soc)
     225             : 
     226        1054 :       IF (do_soc) THEN
     227             :          ! start with multiplicity that is not specified in input
     228             :          ! so that excited-state gradient is for multiplicity given in input
     229           8 :          lmult_tmp = tddfpt_control%rks_triplets
     230           8 :          tddfpt_control%rks_triplets = .NOT. (tddfpt_control%rks_triplets)
     231             :       END IF
     232             : 
     233        1054 :       CALL cite_reference(Iannuzzi2005)
     234        1054 :       IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
     235         402 :          CALL cite_reference(Grimme2013)
     236         402 :          CALL cite_reference(Grimme2016)
     237             :       END IF
     238             : 
     239        1054 :       CALL tddfpt_header(log_unit)
     240        1054 :       CALL kernel_info(log_unit, dft_control, tddfpt_control, xc_section)
     241             :       ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals
     242        1054 :       NULLIFY (gs_mos)
     243             : 
     244        1054 :       CALL tddfpt_init_mos(qs_env, gs_mos, log_unit)
     245             : 
     246             :       ! obtain corrected KS-matrix
     247        1054 :       CALL tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
     248             : 
     249        1054 :       IF ((tddfpt_control%do_lrigpw) .AND. &
     250             :           (tddfpt_control%kernel /= tddfpt_kernel_full)) THEN
     251           0 :          CALL cp_abort(__LOCATION__, "LRI only implemented for full kernel")
     252             :       END IF
     253             : 
     254        1054 :       IF (ASSOCIATED(matrix_ks_oep)) matrix_ks => matrix_ks_oep
     255             : 
     256             :       ! components of the dipole operator
     257             :       CALL tddfpt_dipole_operator(dipole_op_mos_occ, &
     258             :                                   tddfpt_control, &
     259             :                                   gs_mos, &
     260        1054 :                                   qs_env)
     261             : 
     262        1054 :       nspins = SIZE(gs_mos)
     263             :       ! multiplicity of molecular system
     264        1054 :       IF (nspins > 1) THEN
     265         124 :          mult = ABS(SIZE(gs_mos(1)%evals_occ) - SIZE(gs_mos(2)%evals_occ)) + 1
     266         124 :          IF (mult > 2) &
     267           0 :             CALL cp_warn(__LOCATION__, "There is a convergence issue for multiplicity >= 3")
     268             :       ELSE
     269         930 :          IF (tddfpt_control%rks_triplets) THEN
     270         174 :             mult = 3
     271             :          ELSE
     272         756 :             mult = 1
     273             :          END IF
     274             :       END IF
     275             : 
     276             :       ! split mpi communicator
     277        4340 :       ALLOCATE (my_mos(nspins))
     278        2232 :       DO ispin = 1, nspins
     279        2232 :          my_mos(ispin) = gs_mos(ispin)%mos_occ
     280             :       END DO
     281             :       CALL tddfpt_sub_env_init(sub_env, qs_env, mos_occ=my_mos(:), &
     282        1054 :                                kernel=tddfpt_control%kernel)
     283        1054 :       DEALLOCATE (my_mos)
     284             : 
     285        1054 :       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     286             :          ! create environment for Full Kernel
     287         558 :          IF (dft_control%qs_control%xtb) THEN
     288           0 :             CPABORT("TDDFPT: xTB only works with sTDA Kernel")
     289             :          END IF
     290             : 
     291         558 :          IF (tddfpt_control%do_hfxsr) THEN
     292           4 :             kernel_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL")
     293             :             CALL section_vals_val_get(kernel_section, "HFXSR_PRIMBAS", &
     294           4 :                                       i_val=tddfpt_control%hfxsr_primbas)
     295             :             ! basis set
     296             :             CALL create_minbas_set(qs_env, log_unit, basis_type="TDA_HFX", &
     297           4 :                                    primitive=tddfpt_control%hfxsr_primbas)
     298             :             ! admm control
     299          20 :             ALLOCATE (full_kernel_env%admm_control)
     300           4 :             full_kernel_env%admm_control%purification_method = do_admm_purify_none
     301             :             full_kernel_env%admm_control%method = do_admm_basis_projection
     302             :             full_kernel_env%admm_control%scaling_model = do_admm_exch_scaling_none
     303           4 :             full_kernel_env%admm_control%aux_exch_func = do_admm_aux_exch_func_none
     304             :             ! hfx section
     305           4 :             full_kernel_env%hfxsr_section => hfxsr_section
     306             :             !
     307             :             CALL aux_admm_init(qs_env, mos, full_kernel_env%admm_env, &
     308           4 :                                full_kernel_env%admm_control, "TDA_HFX")
     309             :             CALL get_admm_env(full_kernel_env%admm_env, mos_aux_fit=mos_aux_fit, &
     310             :                               matrix_s_aux_fit=matrix_s_aux_fit, &
     311           4 :                               matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
     312             :             CALL admm_fit_mo_coeffs(full_kernel_env%admm_env, matrix_s_aux_fit, &
     313           4 :                                     matrix_s_aux_fit_vs_orb, mos, mos_aux_fit, .TRUE.)
     314             :             ! x_data
     315             :             CALL get_qs_env(qs_env, cell=cell, atomic_kind_set=atomic_kind_set, &
     316             :                             qs_kind_set=qs_kind_set, particle_set=particle_set, &
     317           4 :                             para_env=para_env)
     318             :             CALL hfx_create(full_kernel_env%x_data, para_env, hfxsr_section, atomic_kind_set, &
     319           4 :                             qs_kind_set, particle_set, dft_control, cell, orb_basis="TDA_HFX")
     320             :          END IF
     321             : 
     322             :          ! allocate pools and work matrices
     323         558 :          nstates = tddfpt_control%nstates
     324             :          !! Too many states can lead to Problems
     325             :          !! You should be warned if there are more states
     326             :          !! then occ-virt Combinations!!
     327         558 :          CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
     328         558 :          CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
     329         558 :          nstate_max = nocc*nvirt
     330         558 :          IF (nstates > nstate_max) THEN
     331           0 :             CPWARN("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
     332           0 :             CPWARN("Experimental: CHANGED NSTATES TO ITS MAXIMUM VALUE!")
     333           0 :             nstates = nstate_max
     334           0 :             tddfpt_control%nstates = nstate_max
     335             :          END IF
     336             :          CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, do_admm, &
     337         558 :                                           do_hfxlr, do_exck, qs_env, sub_env)
     338             : 
     339             :          ! create full_kernel and admm_kernel within tddfpt_energies
     340         558 :          kernel_env%full_kernel => full_kernel_env
     341         558 :          kernel_env%admm_kernel => kernel_env_admm_aux
     342         558 :          NULLIFY (kernel_env%stda_kernel)
     343         558 :          IF (do_hfxsr) THEN
     344             :             ! work matrices for SR HFX
     345           4 :             CALL hfxsr_create_work_matrices(work_matrices, qs_env, full_kernel_env%admm_env)
     346             :          END IF
     347         558 :          IF (do_hfxlr) THEN
     348             :             ! calculate S_half and Lowdin MO coefficients
     349           6 :             CALL get_lowdin_mo_coefficients(qs_env, sub_env, work_matrices)
     350             :          END IF
     351         496 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
     352             :          ! setup for kernel_stda outside tddfpt_energies
     353         402 :          nactive = 0
     354         402 :          CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
     355         836 :          DO ispin = 1, SIZE(gs_mos)
     356             :             CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, &
     357         836 :                                 ncol_global=nactive(ispin))
     358             :          END DO
     359         402 :          CALL allocate_stda_env(qs_env, stda_kernel, nao, nactive)
     360             :          ! sTDA parameters
     361         402 :          CALL stda_init_param(qs_env, stda_kernel, tddfpt_control%stda_control)
     362             :          ! allocate pools and work matrices
     363         402 :          nstates = tddfpt_control%nstates
     364             :          CALL stda_create_work_matrices(work_matrices, gs_mos, nstates, &
     365         402 :                                         qs_env, sub_env)
     366             :          !
     367             :          CALL stda_init_matrices(qs_env, stda_kernel, sub_env, &
     368         402 :                                  work_matrices, tddfpt_control)
     369             :          !
     370         402 :          kernel_env%stda_kernel => stda_kernel
     371         402 :          NULLIFY (kernel_env%full_kernel)
     372         402 :          NULLIFY (kernel_env%admm_kernel)
     373          94 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
     374             :          ! allocate pools and work matrices
     375          94 :          nstates = tddfpt_control%nstates
     376             :          CALL stda_create_work_matrices(work_matrices, gs_mos, nstates, &
     377          94 :                                         qs_env, sub_env)
     378          94 :          NULLIFY (kernel_env%full_kernel)
     379          94 :          NULLIFY (kernel_env%admm_kernel)
     380          94 :          NULLIFY (kernel_env%stda_kernel)
     381             :       END IF
     382             : 
     383       10298 :       ALLOCATE (evects(nspins, nstates))
     384        3162 :       ALLOCATE (evals(nstates))
     385             : 
     386       10298 :       ALLOCATE (S_evects(nspins, nstates))
     387        3906 :       DO istate = 1, nstates
     388        7136 :          DO ispin = 1, nspins
     389             :             CALL fm_pool_create_fm( &
     390             :                work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
     391        6082 :                S_evects(ispin, istate))
     392             :          END DO
     393             :       END DO
     394             : 
     395        1054 :       IF (.NOT. do_soc) THEN
     396             :          ! compute tddfpt excitation energies of multiplicity mult
     397             :          CALL tddfpt_energies(qs_env, nstates, work_matrices, &
     398             :                               tddfpt_control, logger, tddfpt_print_section, evects, evals, &
     399             :                               gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
     400             :                               sub_env, ostrength, dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
     401        1046 :                               kernel_env_admm_aux)
     402             :       ELSE
     403             :          CALL tddfpt_soc_energies(qs_env, nstates, work_matrices, &
     404             :                                   tddfpt_control, logger, tddfpt_print_section, &
     405             :                                   evects, evals, ostrength, &
     406             :                                   gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
     407             :                                   sub_env, dipole_op_mos_occ, lmult_tmp, xc_section, full_kernel_env, &
     408           8 :                                   kernel_env_admm_aux)
     409             :       END IF
     410             : 
     411             :       !print forces for selected states
     412        1054 :       IF (calc_forces) THEN
     413             :          CALL tddfpt_print_forces(qs_env, evects, evals, ostrength, &
     414             :                                   tddfpt_print_section, gs_mos, &
     415         552 :                                   kernel_env, sub_env, work_matrices)
     416             :       END IF
     417             : 
     418             :       ! excited state potential energy surface
     419        1054 :       IF (qs_env%excited_state) THEN
     420         916 :          IF (sub_env%is_split) THEN
     421             :             CALL cp_abort(__LOCATION__, &
     422             :                           "Excited state forces not possible when states"// &
     423           0 :                           " are distributed to different CPU pools.")
     424             :          END IF
     425             :          ! for gradients unshifted KS matrix
     426         916 :          IF (ASSOCIATED(matrix_ks_oep)) CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
     427         916 :          CALL get_qs_env(qs_env, exstate_env=ex_env)
     428         916 :          state_change = .FALSE.
     429         916 :          IF (ex_env%state > 0) THEN
     430         908 :             my_state = ex_env%state
     431           8 :          ELSEIF (ex_env%state < 0) THEN
     432             :             ! state following
     433          32 :             ALLOCATE (my_mos(nspins))
     434          16 :             DO ispin = 1, nspins
     435          16 :                my_mos(ispin) = gs_mos(ispin)%mos_occ
     436             :             END DO
     437           8 :             my_state = ABS(ex_env%state)
     438           8 :             CALL assign_state(qs_env, matrix_s, evects, my_mos, ex_env%wfn_history, my_state)
     439           8 :             DEALLOCATE (my_mos)
     440           8 :             IF (my_state /= ABS(ex_env%state)) THEN
     441           0 :                state_change = .TRUE.
     442           0 :                old_state = ABS(ex_env%state)
     443             :             END IF
     444           8 :             ex_env%state = -my_state
     445             :          ELSE
     446             :             CALL cp_warn(__LOCATION__, &
     447           0 :                          "Active excited state not assigned. Use the first state.")
     448           0 :             my_state = 1
     449             :          END IF
     450         916 :          CPASSERT(my_state > 0)
     451         916 :          IF (my_state > nstates) THEN
     452             :             CALL cp_warn(__LOCATION__, &
     453           0 :                          "There were not enough excited states calculated.")
     454           0 :             CPABORT("excited state potential energy surface")
     455             :          END IF
     456             :          !
     457             :          ! energy
     458         916 :          ex_env%evalue = evals(my_state)
     459             :          ! excitation vector
     460         916 :          CALL cp_fm_release(ex_env%evect)
     461        3760 :          ALLOCATE (ex_env%evect(nspins))
     462        1928 :          DO ispin = 1, nspins
     463             :             CALL cp_fm_get_info(matrix=evects(ispin, 1), &
     464        1012 :                                 matrix_struct=matrix_struct)
     465        1012 :             CALL cp_fm_create(ex_env%evect(ispin), matrix_struct)
     466        1928 :             CALL cp_fm_to_fm(evects(ispin, my_state), ex_env%evect(ispin))
     467             :          END DO
     468             : 
     469         916 :          IF (log_unit > 0) THEN
     470         458 :             gsval = ex_env%wfn_history%gsval
     471         458 :             gsmin = ex_env%wfn_history%gsmin
     472         458 :             xsval = ex_env%wfn_history%xsval
     473         458 :             WRITE (log_unit, "(1X,A,T40,F10.6,A,T62,F10.6,A)") "Ground state orbital alignment:", &
     474         916 :                gsmin, "[MinVal]", gsval, "[Average]"
     475         458 :             WRITE (log_unit, "(1X,A,T71,F10.6)") "Excitation vector alignment:", xsval
     476         458 :             IF (state_change) THEN
     477             :                WRITE (log_unit, "(1X,A,I5,T60,A14,T76,I5)") &
     478           0 :                   "Target state has been changed from state ", &
     479           0 :                   old_state, " to new state ", my_state
     480             :             END IF
     481         458 :             WRITE (log_unit, "(1X,A,I4,A,F12.5,A)") "Calculate properties for state:", &
     482         916 :                my_state, "      with excitation energy ", ex_env%evalue*evolt, " eV"
     483             :          END IF
     484             : 
     485         916 :          IF (calc_forces) THEN
     486             :             CALL tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, &
     487         550 :                                     sub_env, work_matrices)
     488             :          END IF
     489             :       END IF
     490             : 
     491             :       ! clean up
     492        1054 :       CALL cp_fm_release(evects)
     493        1054 :       CALL cp_fm_release(S_evects)
     494             : 
     495             :       CALL cp_print_key_finished_output(log_unit, &
     496             :                                         logger, &
     497             :                                         tddfpt_print_section, &
     498        1054 :                                         "PROGRAM_BANNER")
     499             : 
     500        1054 :       DEALLOCATE (evals, ostrength)
     501             : 
     502        1054 :       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     503         558 :          IF (do_admm) CALL release_kernel_env(kernel_env%admm_kernel)
     504         558 :          IF (tddfpt_control%do_lrigpw) THEN
     505          10 :             CALL lri_env_release(kernel_env%full_kernel%lri_env)
     506          10 :             DEALLOCATE (kernel_env%full_kernel%lri_env)
     507          10 :             CALL lri_density_release(kernel_env%full_kernel%lri_density)
     508          10 :             DEALLOCATE (kernel_env%full_kernel%lri_density)
     509             :          END IF
     510         558 :          CALL release_kernel_env(kernel_env%full_kernel)
     511         496 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
     512         402 :          CALL deallocate_stda_env(stda_kernel)
     513          94 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
     514             :          !
     515             :       ELSE
     516           0 :          CPABORT('Unknown kernel type')
     517             :       END IF
     518        1054 :       CALL tddfpt_release_work_matrices(work_matrices, sub_env)
     519        1054 :       CALL tddfpt_sub_env_release(sub_env)
     520             : 
     521        1054 :       CALL cp_fm_release(dipole_op_mos_occ)
     522             : 
     523        2232 :       DO ispin = nspins, 1, -1
     524        2232 :          CALL tddfpt_release_ground_state_mos(gs_mos(ispin))
     525             :       END DO
     526        1054 :       DEALLOCATE (gs_mos)
     527             : 
     528        1054 :       IF (ASSOCIATED(matrix_ks_oep)) &
     529          30 :          CALL dbcsr_deallocate_matrix_set(matrix_ks_oep)
     530             : 
     531        1054 :       CALL timestop(handle)
     532             : 
     533        8432 :    END SUBROUTINE tddfpt
     534             : 
     535             : ! **************************************************************************************************
     536             : !> \brief TDDFPT input
     537             : !> \param qs_env  Quickstep environment
     538             : !> \param do_hfx ...
     539             : !> \param do_admm ...
     540             : !> \param do_exck ...
     541             : !> \param do_hfxsr ...
     542             : !> \param do_hfxlr ...
     543             : !> \param xc_section ...
     544             : !> \param tddfpt_print_section ...
     545             : !> \param lri_section ...
     546             : !> \param hfxsr_section ...
     547             : ! **************************************************************************************************
     548        1054 :    SUBROUTINE tddfpt_input(qs_env, do_hfx, do_admm, do_exck, do_hfxsr, do_hfxlr, &
     549             :                            xc_section, tddfpt_print_section, lri_section, hfxsr_section)
     550             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     551             :       LOGICAL, INTENT(INOUT)                             :: do_hfx, do_admm, do_exck, do_hfxsr, &
     552             :                                                             do_hfxlr
     553             :       TYPE(section_vals_type), POINTER                   :: xc_section, tddfpt_print_section, &
     554             :                                                             lri_section, hfxsr_section
     555             : 
     556             :       CHARACTER(len=20)                                  :: nstates_str
     557             :       LOGICAL                                            :: exar, exf, exgcp, exhf, exhfxk, exk, &
     558             :                                                             explicit_root, expot, exvdw, exwfn, &
     559             :                                                             found, same_hfx
     560             :       REAL(kind=dp)                                      :: C_hf
     561             :       TYPE(dft_control_type), POINTER                    :: dft_control
     562             :       TYPE(section_vals_type), POINTER                   :: hfx_section, hfx_section_gs, input, &
     563             :                                                             tddfpt_section, xc_root, xc_sub
     564             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     565             : 
     566        1054 :       NULLIFY (dft_control, input)
     567        1054 :       CALL get_qs_env(qs_env, dft_control=dft_control, input=input)
     568        1054 :       tddfpt_control => dft_control%tddfpt2_control
     569             : 
     570             :       ! no k-points
     571        1054 :       CPASSERT(dft_control%nimages <= 1)
     572             : 
     573        1054 :       IF (tddfpt_control%nstates <= 0) THEN
     574           0 :          CALL integer_to_string(tddfpt_control%nstates, nstates_str)
     575             :          CALL cp_warn(__LOCATION__, "TDDFPT calculation was requested for "// &
     576           0 :                       TRIM(nstates_str)//" excited states: nothing to do.")
     577           0 :          RETURN
     578             :       END IF
     579             : 
     580        1054 :       NULLIFY (tddfpt_section, tddfpt_print_section)
     581        1054 :       tddfpt_section => section_vals_get_subs_vals(input, "PROPERTIES%TDDFPT")
     582        1054 :       tddfpt_print_section => section_vals_get_subs_vals(tddfpt_section, "PRINT")
     583             : 
     584        1054 :       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     585         558 :          NULLIFY (xc_root)
     586         558 :          xc_root => section_vals_get_subs_vals(tddfpt_section, "XC")
     587         558 :          CALL section_vals_get(xc_root, explicit=explicit_root)
     588         558 :          NULLIFY (xc_section)
     589         558 :          IF (explicit_root) THEN
     590             :             ! No ADIABATIC_RESCALING option possible
     591         334 :             NULLIFY (xc_sub)
     592         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "ADIABATIC_RESCALING")
     593         334 :             CALL section_vals_get(xc_sub, explicit=exar)
     594         334 :             IF (exar) THEN
     595           0 :                CALL cp_warn(__LOCATION__, "TDDFPT Kernel with ADIABATIC_RESCALING not possible.")
     596           0 :                CPABORT("TDDFPT Input")
     597             :             END IF
     598             :             ! No GCP_POTENTIAL option possible
     599         334 :             NULLIFY (xc_sub)
     600         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "GCP_POTENTIAL")
     601         334 :             CALL section_vals_get(xc_sub, explicit=exgcp)
     602         334 :             IF (exgcp) THEN
     603           0 :                CALL cp_warn(__LOCATION__, "TDDFPT Kernel with GCP_POTENTIAL not possible.")
     604           0 :                CPABORT("TDDFPT Input")
     605             :             END IF
     606             :             ! No VDW_POTENTIAL option possible
     607         334 :             NULLIFY (xc_sub)
     608         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "VDW_POTENTIAL")
     609         334 :             CALL section_vals_get(xc_sub, explicit=exvdw)
     610         334 :             IF (exvdw) THEN
     611           0 :                CALL cp_warn(__LOCATION__, "TDDFPT Kernel with VDW_POTENTIAL not possible.")
     612           0 :                CPABORT("TDDFPT Input")
     613             :             END IF
     614             :             ! No WF_CORRELATION option possible
     615         334 :             NULLIFY (xc_sub)
     616         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "WF_CORRELATION")
     617         334 :             CALL section_vals_get(xc_sub, explicit=exwfn)
     618         334 :             IF (exwfn) THEN
     619           0 :                CALL cp_warn(__LOCATION__, "TDDFPT Kernel with WF_CORRELATION not possible.")
     620           0 :                CPABORT("TDDFPT Input")
     621             :             END IF
     622             :             ! No XC_POTENTIAL option possible
     623         334 :             NULLIFY (xc_sub)
     624         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "XC_POTENTIAL")
     625         334 :             CALL section_vals_get(xc_sub, explicit=expot)
     626         334 :             IF (expot) THEN
     627           0 :                CALL cp_warn(__LOCATION__, "TDDFPT Kernel with XC_POTENTIAL not possible.")
     628           0 :                CPABORT("TDDFPT Input")
     629             :             END IF
     630             :             !
     631         334 :             NULLIFY (xc_sub)
     632         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "XC_FUNCTIONAL")
     633         334 :             CALL section_vals_get(xc_sub, explicit=exf)
     634         334 :             NULLIFY (xc_sub)
     635         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "XC_KERNEL")
     636         334 :             CALL section_vals_get(xc_sub, explicit=exk)
     637         334 :             IF ((exf .AND. exk) .OR. .NOT. (exf .OR. exk)) THEN
     638           0 :                CALL cp_warn(__LOCATION__, "TDDFPT Kernel needs XC_FUNCTIONAL or XC_KERNEL section.")
     639           0 :                CPABORT("TDDFPT Input")
     640             :             END IF
     641         334 :             NULLIFY (xc_sub)
     642         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "HF")
     643         334 :             CALL section_vals_get(xc_sub, explicit=exhf)
     644         334 :             NULLIFY (xc_sub)
     645         334 :             xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL")
     646         334 :             CALL section_vals_get(xc_sub, explicit=exhfxk)
     647             :             !
     648         334 :             xc_section => xc_root
     649         334 :             hfx_section => section_vals_get_subs_vals(xc_section, "HF")
     650         334 :             CALL section_vals_get(hfx_section, explicit=do_hfx)
     651         334 :             IF (do_hfx) THEN
     652          12 :                CALL section_vals_val_get(hfx_section, "FRACTION", r_val=C_hf)
     653          12 :                do_hfx = (C_hf /= 0.0_dp)
     654             :             END IF
     655             :             !TDDFPT only works if the kernel has the same HF section as the DFT%XC one
     656         334 :             IF (do_hfx) THEN
     657          12 :                hfx_section_gs => section_vals_get_subs_vals(input, "DFT%XC%HF")
     658          12 :                CALL compare_hfx_sections(hfx_section, hfx_section_gs, same_hfx)
     659          12 :                IF (.NOT. same_hfx) THEN
     660           0 :                   CPABORT("TDDFPT Kernel must use the same HF section as DFT%XC or no HF at all.")
     661             :                END IF
     662             :             END IF
     663             : 
     664         334 :             do_admm = do_hfx .AND. dft_control%do_admm
     665         334 :             IF (do_admm) THEN
     666             :                ! 'admm_env%xc_section_primary' and 'admm_env%xc_section_aux' need to be redefined
     667             :                CALL cp_abort(__LOCATION__, &
     668             :                              "ADMM is not implemented for a TDDFT kernel XC-functional which is different from "// &
     669           0 :                              "the one used for the ground-state calculation. A ground-state 'admm_env' cannot be reused.")
     670             :             END IF
     671             :             ! SET HFX_KERNEL and/or XC_KERNEL
     672         334 :             IF (exk) THEN
     673          12 :                do_exck = .TRUE.
     674             :             ELSE
     675         322 :                do_exck = .FALSE.
     676             :             END IF
     677         334 :             IF (exhfxk) THEN
     678           6 :                xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL")
     679           6 :                CALL section_vals_val_get(xc_sub, "DO_HFXSR", l_val=do_hfxsr)
     680           6 :                xc_sub => section_vals_get_subs_vals(xc_root, "HFX_KERNEL%HFXLR")
     681           6 :                CALL section_vals_get(xc_sub, explicit=do_hfxlr)
     682             :             ELSE
     683         328 :                do_hfxsr = .FALSE.
     684         328 :                do_hfxlr = .FALSE.
     685             :             END IF
     686             :          ELSE
     687         224 :             xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     688         224 :             hfx_section => section_vals_get_subs_vals(xc_section, "HF")
     689         224 :             CALL section_vals_get(hfx_section, explicit=do_hfx)
     690         224 :             IF (do_hfx) THEN
     691         194 :                CALL section_vals_val_get(hfx_section, "FRACTION", r_val=C_hf)
     692         194 :                do_hfx = (C_hf /= 0.0_dp)
     693             :             END IF
     694         224 :             do_admm = do_hfx .AND. dft_control%do_admm
     695         224 :             do_exck = .FALSE.
     696         224 :             do_hfxsr = .FALSE.
     697         224 :             do_hfxlr = .FALSE.
     698             :          END IF
     699             :       ELSE
     700         496 :          do_hfx = .FALSE.
     701         496 :          do_admm = .FALSE.
     702         496 :          do_exck = .FALSE.
     703         496 :          do_hfxsr = .FALSE.
     704         496 :          do_hfxlr = .FALSE.
     705             :       END IF
     706             : 
     707             :       ! reset rks_triplets if UKS is in use
     708        1054 :       IF (tddfpt_control%rks_triplets .AND. dft_control%nspins > 1) THEN
     709           8 :          tddfpt_control%rks_triplets = .FALSE.
     710           8 :          CALL cp_warn(__LOCATION__, "Keyword RKS_TRIPLETS has been ignored for spin-polarised calculations")
     711             :       END IF
     712             : 
     713             :       ! lri input
     714        1054 :       IF (tddfpt_control%do_lrigpw) THEN
     715          10 :          lri_section => section_vals_get_subs_vals(tddfpt_section, "LRIGPW")
     716             :       END IF
     717             : 
     718             :       ! set defaults for short range HFX
     719        1054 :       NULLIFY (hfxsr_section)
     720        1054 :       IF (do_hfxsr) THEN
     721           4 :          hfxsr_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL%HF")
     722           4 :          CALL section_vals_get(hfxsr_section, explicit=found)
     723           4 :          IF (.NOT. found) THEN
     724           0 :             CPABORT("HFXSR option needs &HF section defined")
     725             :          END IF
     726           4 :          CALL section_vals_val_get(hfxsr_section, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", explicit=found)
     727           4 :          IF (.NOT. found) THEN
     728             :             CALL section_vals_val_set(hfxsr_section, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", &
     729           4 :                                       i_val=do_potential_truncated)
     730             :          END IF
     731           4 :          CALL section_vals_val_get(hfxsr_section, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", explicit=found)
     732           4 :          IF (.NOT. found) THEN
     733           4 :             CALL section_vals_val_set(hfxsr_section, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=7.5589_dp)
     734             :          END IF
     735           4 :          CALL section_vals_val_get(hfxsr_section, "RI%_SECTION_PARAMETERS_", l_val=found)
     736           4 :          IF (found) THEN
     737           0 :             CALL cp_abort(__LOCATION__, "Short range TDA kernel with RI not possible")
     738             :          END IF
     739             :       END IF
     740             : 
     741             :    END SUBROUTINE tddfpt_input
     742             : 
     743             : ! **************************************************************************************************
     744             : !> \brief ...
     745             : !> \param log_unit ...
     746             : !> \param dft_control ...
     747             : !> \param tddfpt_control ...
     748             : !> \param xc_section ...
     749             : ! **************************************************************************************************
     750        1054 :    SUBROUTINE kernel_info(log_unit, dft_control, tddfpt_control, xc_section)
     751             :       INTEGER, INTENT(IN)                                :: log_unit
     752             :       TYPE(dft_control_type), POINTER                    :: dft_control
     753             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     754             :       TYPE(section_vals_type), POINTER                   :: xc_section
     755             : 
     756             :       CHARACTER(LEN=4)                                   :: ktype
     757             :       LOGICAL                                            :: lsd
     758             : 
     759        1054 :       lsd = (dft_control%nspins > 1)
     760        1054 :       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     761         558 :          ktype = "FULL"
     762         558 :          IF (log_unit > 0) THEN
     763         279 :             WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", TRIM(ktype)
     764         279 :             CALL xc_write(log_unit, xc_section, lsd)
     765         279 :             IF (tddfpt_control%do_hfx) THEN
     766         103 :                IF (tddfpt_control%do_admm) THEN
     767          60 :                   WRITE (log_unit, "(T2,A,T62,A19)") "KERNEL|", "ADMM Exact Exchange"
     768          60 :                   IF (tddfpt_control%admm_xc_correction) THEN
     769          46 :                      WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Apply ADMM Kernel XC Correction"
     770             :                   END IF
     771          60 :                   IF (tddfpt_control%admm_symm) THEN
     772          60 :                      WRITE (log_unit, "(T2,A,T60,A21)") "KERNEL|", "Symmetric ADMM Kernel"
     773             :                   END IF
     774             :                ELSE
     775          43 :                   WRITE (log_unit, "(T2,A,T67,A14)") "KERNEL|", "Exact Exchange"
     776             :                END IF
     777             :             END IF
     778         279 :             IF (tddfpt_control%do_hfxsr) THEN
     779           2 :                WRITE (log_unit, "(T2,A,T43,A38)") "KERNEL|", "Short range HFX approximation"
     780             :             END IF
     781         279 :             IF (tddfpt_control%do_hfxlr) THEN
     782           3 :                WRITE (log_unit, "(T2,A,T43,A38)") "KERNEL|", "Long range HFX approximation"
     783             :             END IF
     784         279 :             IF (tddfpt_control%do_lrigpw) THEN
     785           5 :                WRITE (log_unit, "(T2,A,T42,A39)") "KERNEL|", "LRI approximation of transition density"
     786             :             END IF
     787             :          END IF
     788         496 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
     789         402 :          ktype = "sTDA"
     790         402 :          IF (log_unit > 0) THEN
     791         201 :             WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", TRIM(ktype)
     792         201 :             IF (tddfpt_control%stda_control%do_ewald) THEN
     793          47 :                WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Coulomb term uses Ewald summation"
     794             :             ELSE
     795         154 :                WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Coulomb term uses direct summation (MIC)"
     796             :             END IF
     797         201 :             IF (tddfpt_control%stda_control%do_exchange) THEN
     798         185 :                WRITE (log_unit, "(T2,A,T78,A3)") "KERNEL| Exact exchange term", "YES"
     799         185 :                WRITE (log_unit, "(T2,A,T71,F10.3)") "KERNEL| Short range HFX fraction:", &
     800         370 :                   tddfpt_control%stda_control%hfx_fraction
     801             :             ELSE
     802          16 :                WRITE (log_unit, "(T2,A,T79,A2)") "KERNEL| Exact exchange term", "NO"
     803             :             END IF
     804         201 :             WRITE (log_unit, "(T2,A,T66,E15.3)") "KERNEL| Transition density filter", &
     805         402 :                tddfpt_control%stda_control%eps_td_filter
     806             :          END IF
     807          94 :       ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
     808          94 :          ktype = "NONE"
     809          94 :          IF (log_unit > 0) THEN
     810          47 :             WRITE (log_unit, "(T2,A,T77,A4)") "KERNEL|", TRIM(ktype)
     811             :          END IF
     812             :       ELSE
     813             :          !CPABORT("Unknown kernel")
     814             :       END IF
     815             :       !
     816        1054 :       IF (log_unit > 0) THEN
     817         527 :          IF (tddfpt_control%rks_triplets) THEN
     818          87 :             WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Triplet"
     819         440 :          ELSE IF (lsd) THEN
     820          62 :             WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin symmetry of excitations", "Unrestricted"
     821             :          ELSE
     822         378 :             WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Singlet"
     823             :          END IF
     824         527 :          WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Number of states calculated", tddfpt_control%nstates
     825         527 :          WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Number of Davidson iterations", tddfpt_control%niters
     826         527 :          WRITE (log_unit, "(T2,A,T66,E15.3)") "TDDFPT| Davidson iteration convergence", tddfpt_control%conv
     827         527 :          WRITE (log_unit, "(T2,A,T73,I8)") "TDDFPT| Max. number of Krylov space vectors", tddfpt_control%nkvs
     828             :       END IF
     829             : 
     830        1054 :    END SUBROUTINE kernel_info
     831             : 
     832             : ! **************************************************************************************************
     833             : !> \brief The energy cylculation has been moved to its own subroutine
     834             : !> \param qs_env ...
     835             : !> \param nstates ...
     836             : !> \param work_matrices ...
     837             : !> \param tddfpt_control ...
     838             : !> \param logger ...
     839             : !> \param tddfpt_print_section ...
     840             : !> \param evects ...
     841             : !> \param evals ...
     842             : !> \param gs_mos ...
     843             : !> \param tddfpt_section ...
     844             : !> \param S_evects ...
     845             : !> \param matrix_s ...
     846             : !> \param kernel_env ...
     847             : !> \param matrix_ks ...
     848             : !> \param sub_env ...
     849             : !> \param ostrength ...
     850             : !> \param dipole_op_mos_occ ...
     851             : !> \param mult ...
     852             : !> \param xc_section ...
     853             : !> \param full_kernel_env ...
     854             : !> \param kernel_env_admm_aux ...
     855             : ! **************************************************************************************************
     856        1062 :    SUBROUTINE tddfpt_energies(qs_env, nstates, work_matrices, &
     857             :                               tddfpt_control, logger, tddfpt_print_section, evects, evals, &
     858             :                               gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
     859             :                               sub_env, ostrength, dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
     860             :                               kernel_env_admm_aux)
     861             : 
     862             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     863             :       INTEGER                                            :: nstates
     864             :       TYPE(tddfpt_work_matrices)                         :: work_matrices
     865             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     866             :       TYPE(cp_logger_type), POINTER                      :: logger
     867             :       TYPE(section_vals_type), POINTER                   :: tddfpt_print_section
     868             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: evects
     869             :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: evals
     870             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     871             :          POINTER                                         :: gs_mos
     872             :       TYPE(section_vals_type), POINTER                   :: tddfpt_section
     873             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: S_evects
     874             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     875             :       TYPE(kernel_env_type)                              :: kernel_env
     876             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     877             :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
     878             :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: ostrength
     879             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: dipole_op_mos_occ
     880             :       INTEGER                                            :: mult
     881             :       TYPE(section_vals_type), POINTER                   :: xc_section
     882             :       TYPE(full_kernel_env_type), TARGET                 :: full_kernel_env, kernel_env_admm_aux
     883             : 
     884             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_energies'
     885             : 
     886             :       CHARACTER(len=20)                                  :: nstates_str
     887             :       INTEGER                                            :: energy_unit, handle, iter, log_unit, &
     888             :                                                             niters, nocc, nstate_max, &
     889             :                                                             nstates_read, nvirt
     890             :       LOGICAL                                            :: do_admm, do_exck, do_soc, explicit, &
     891             :                                                             is_restarted
     892             :       REAL(kind=dp)                                      :: conv
     893             :       TYPE(admm_type), POINTER                           :: admm_env
     894             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     895        1062 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks_oep
     896             :       TYPE(section_vals_type), POINTER                   :: lri_section, namd_print_section, &
     897             :                                                             soc_section
     898             : 
     899        1062 :       CALL timeset(routineN, handle)
     900             : 
     901        1062 :       NULLIFY (admm_env, matrix_ks_oep)
     902        1062 :       do_admm = tddfpt_control%do_admm
     903        1062 :       IF (do_admm) CALL get_qs_env(qs_env, admm_env=admm_env)
     904             : 
     905             :       ! setup for full_kernel and admm_kernel within tddfpt_energies due to dependence on multiplicity
     906        1062 :       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
     907             : 
     908             :          CALL tddfpt_construct_ground_state_orb_density( &
     909             :             rho_orb_struct=work_matrices%rho_orb_struct_sub, &
     910             :             rho_xc_struct=work_matrices%rho_xc_struct_sub, &
     911             :             is_rks_triplets=tddfpt_control%rks_triplets, &
     912             :             qs_env=qs_env, sub_env=sub_env, &
     913         566 :             wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub)
     914             : 
     915         566 :          IF (do_admm) THEN
     916             :             ! Full kernel with ADMM
     917         120 :             IF (tddfpt_control%admm_xc_correction) THEN
     918             :                CALL create_kernel_env(kernel_env=full_kernel_env, &
     919             :                                       rho_struct_sub=work_matrices%rho_orb_struct_sub, &
     920             :                                       xc_section=admm_env%xc_section_primary, &
     921             :                                       is_rks_triplets=tddfpt_control%rks_triplets, &
     922          92 :                                       sub_env=sub_env)
     923             :             ELSE
     924             :                CALL create_kernel_env(kernel_env=full_kernel_env, &
     925             :                                       rho_struct_sub=work_matrices%rho_orb_struct_sub, &
     926             :                                       xc_section=xc_section, &
     927             :                                       is_rks_triplets=tddfpt_control%rks_triplets, &
     928          28 :                                       sub_env=sub_env)
     929             :             END IF
     930             : 
     931             :             CALL tddfpt_construct_aux_fit_density( &
     932             :                rho_orb_struct=work_matrices%rho_orb_struct_sub, &
     933             :                rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, &
     934             :                local_rho_set=sub_env%local_rho_set_admm, &
     935             :                qs_env=qs_env, sub_env=sub_env, &
     936             :                wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub, &
     937             :                wfm_rho_aux_fit=work_matrices%rho_ao_aux_fit_fm_sub, &
     938         120 :                wfm_aux_orb=work_matrices%wfm_aux_orb_sub)
     939             : 
     940             :             CALL create_kernel_env(kernel_env=kernel_env_admm_aux, &
     941             :                                    rho_struct_sub=work_matrices%rho_aux_fit_struct_sub, &
     942             :                                    xc_section=admm_env%xc_section_aux, &
     943             :                                    is_rks_triplets=tddfpt_control%rks_triplets, &
     944         120 :                                    sub_env=sub_env)
     945         120 :             kernel_env%full_kernel => full_kernel_env
     946         120 :             kernel_env%admm_kernel => kernel_env_admm_aux
     947             :          ELSE
     948             :             ! Full kernel
     949             :             CALL create_kernel_env(kernel_env=full_kernel_env, &
     950             :                                    rho_struct_sub=work_matrices%rho_orb_struct_sub, &
     951             :                                    xc_section=xc_section, &
     952             :                                    is_rks_triplets=tddfpt_control%rks_triplets, &
     953         446 :                                    sub_env=sub_env)
     954         446 :             kernel_env%full_kernel => full_kernel_env
     955         446 :             NULLIFY (kernel_env%admm_kernel)
     956             :          END IF
     957             :          ! Fxc from kernel definition
     958         566 :          do_exck = tddfpt_control%do_exck
     959         566 :          kernel_env%full_kernel%do_exck = do_exck
     960             :          ! initilize xc kernel
     961         566 :          IF (do_exck) THEN
     962             :             CALL create_fxc_kernel(work_matrices%rho_orb_struct_sub, work_matrices%fxc_rspace_sub, &
     963          12 :                                    xc_section, tddfpt_control%rks_triplets, sub_env, qs_env)
     964             :          END IF
     965             :       END IF
     966             : 
     967             :       ! lri input
     968        1062 :       IF (tddfpt_control%do_lrigpw) THEN
     969          10 :          lri_section => section_vals_get_subs_vals(tddfpt_section, "LRIGPW")
     970             :          CALL tddfpt2_lri_init(qs_env, kernel_env, lri_section, &
     971          10 :                                tddfpt_print_section)
     972             :       END IF
     973             : 
     974             :       !! Too many states cal lead to Problems
     975             :       !! You should be warned if there are mor states
     976             :       !! then occ-virt Combinations!!
     977        1062 :       CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
     978        1062 :       CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
     979        1062 :       nstate_max = nocc*nvirt
     980        1062 :       IF (nstates > nstate_max) THEN
     981           0 :          CPWARN("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
     982           0 :          CPWARN("Experimental: CHANGED NSTATES TO ITS MAXIMUM VALUE!")
     983           0 :          nstates = nstate_max
     984             :       END IF
     985             : 
     986        1062 :       soc_section => section_vals_get_subs_vals(tddfpt_section, "SOC")
     987        1062 :       CALL section_vals_get(soc_section, explicit=do_soc)
     988             : 
     989             :       ! reuse Ritz vectors from the previous calculation if available
     990        1062 :       IF (tddfpt_control%is_restart .AND. .NOT. do_soc) THEN
     991           2 :          CALL get_qs_env(qs_env, blacs_env=blacs_env)
     992             : 
     993             :          nstates_read = tddfpt_read_restart( &
     994             :                         evects=evects, &
     995             :                         evals=evals, &
     996             :                         gs_mos=gs_mos, &
     997             :                         logger=logger, &
     998             :                         tddfpt_section=tddfpt_section, &
     999             :                         tddfpt_print_section=tddfpt_print_section, &
    1000             :                         fm_pool_ao_mo_occ=work_matrices%fm_pool_ao_mo_occ, &
    1001           2 :                         blacs_env_global=blacs_env)
    1002             :       ELSE
    1003             :          nstates_read = 0
    1004             :       END IF
    1005             : 
    1006        1062 :       is_restarted = nstates_read >= nstates
    1007             : 
    1008             :       ! build the list of missed singly excited states and sort them in ascending order
    1009             :       ! according to their excitation energies
    1010             :       log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
    1011        1062 :                                       "GUESS_VECTORS", extension=".tddfptLog")
    1012             :       CALL tddfpt_guess_vectors(evects=evects, evals=evals, &
    1013        1062 :                                 gs_mos=gs_mos, log_unit=log_unit)
    1014             :       CALL cp_print_key_finished_output(log_unit, logger, &
    1015        1062 :                                         tddfpt_print_section, "GUESS_VECTORS")
    1016             : 
    1017        1062 :       CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T)
    1018             :       CALL tddfpt_orthonormalize_psi1_psi1(evects, &
    1019             :                                            nstates, &
    1020             :                                            S_evects, &
    1021        1062 :                                            matrix_s(1)%matrix)
    1022             : 
    1023        1062 :       niters = tddfpt_control%niters
    1024        1062 :       IF (niters > 0) THEN
    1025             :          log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
    1026        1062 :                                          "ITERATION_INFO", extension=".tddfptLog")
    1027             :          energy_unit = cp_print_key_unit_nr(logger, &
    1028             :                                             tddfpt_print_section, &
    1029             :                                             "DETAILED_ENERGY", &
    1030        1062 :                                             extension=".tddfptLog")
    1031             : 
    1032        1062 :          IF (log_unit > 0) THEN
    1033         531 :             WRITE (log_unit, "(1X,A)") "", &
    1034         531 :                "-------------------------------------------------------------------------------", &
    1035         531 :                "-                      TDDFPT WAVEFUNCTION OPTIMIZATION                       -", &
    1036        1062 :                "-------------------------------------------------------------------------------"
    1037             : 
    1038         531 :             WRITE (log_unit, '(/,T11,A,T27,A,T40,A,T62,A)') "Step", "Time", "Convergence", "Conv. states"
    1039         531 :             WRITE (log_unit, '(1X,79("-"))')
    1040             :          END IF
    1041             : 
    1042        1062 :          CALL cp_add_iter_level(logger%iter_info, "TDDFT_SCF")
    1043             : 
    1044             :          DO
    1045             :             ! *** perform Davidson iterations ***
    1046             :             conv = tddfpt_davidson_solver( &
    1047             :                    evects=evects, &
    1048             :                    evals=evals, &
    1049             :                    S_evects=S_evects, &
    1050             :                    gs_mos=gs_mos, &
    1051             :                    tddfpt_control=tddfpt_control, &
    1052             :                    matrix_ks=matrix_ks, &
    1053             :                    qs_env=qs_env, &
    1054             :                    kernel_env=kernel_env, &
    1055             :                    sub_env=sub_env, &
    1056             :                    logger=logger, &
    1057             :                    iter_unit=log_unit, &
    1058             :                    energy_unit=energy_unit, &
    1059             :                    tddfpt_print_section=tddfpt_print_section, &
    1060        2122 :                    work_matrices=work_matrices)
    1061             : 
    1062             :             ! at this point at least one of the following conditions are met:
    1063             :             ! a) convergence criteria has been achieved;
    1064             :             ! b) maximum number of iterations has been reached;
    1065             :             ! c) Davidson iterations must be restarted due to lack of Krylov vectors or numerical instability
    1066             : 
    1067        2122 :             CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter)
    1068             :             ! terminate the loop if either (a) or (b) is true ...
    1069             :             IF ((conv <= tddfpt_control%conv &
    1070        2122 :                  .AND. is_restarted) .OR. iter >= niters) EXIT
    1071             : 
    1072             :             ! ... otherwise restart Davidson iterations
    1073        1060 :             is_restarted = .TRUE.
    1074        2122 :             IF (log_unit > 0) THEN
    1075         530 :                WRITE (log_unit, '(1X,25("-"),1X,A,1X,25("-"))') "Restart Davidson iterations"
    1076         530 :                CALL m_flush(log_unit)
    1077             :             END IF
    1078             :          END DO
    1079             : 
    1080             :          ! write TDDFPT restart file at the last iteration if requested to do so
    1081        1062 :          CALL cp_iterate(logger%iter_info, increment=0, last=.TRUE.)
    1082             :          CALL tddfpt_write_restart(evects=evects, &
    1083             :                                    evals=evals, &
    1084             :                                    gs_mos=gs_mos, &
    1085             :                                    logger=logger, &
    1086        1062 :                                    tddfpt_print_section=tddfpt_print_section)
    1087             : 
    1088        1062 :          CALL cp_rm_iter_level(logger%iter_info, "TDDFT_SCF")
    1089             : 
    1090             :          ! print convergence summary
    1091        1062 :          IF (log_unit > 0) THEN
    1092         531 :             CALL integer_to_string(iter, nstates_str)
    1093         531 :             IF (conv <= tddfpt_control%conv) THEN
    1094         531 :                WRITE (log_unit, "(1X,A)") "", &
    1095         531 :                   "-------------------------------------------------------------------------------", &
    1096         531 :                   "-  TDDFPT run converged in "//TRIM(nstates_str)//" iteration(s) ", &
    1097        1062 :                   "-------------------------------------------------------------------------------"
    1098             :             ELSE
    1099           0 :                WRITE (log_unit, "(1X,A)") "", &
    1100           0 :                   "-------------------------------------------------------------------------------", &
    1101           0 :                   "-  TDDFPT run did NOT converge after "//TRIM(nstates_str)//" iteration(s) ", &
    1102           0 :                   "-------------------------------------------------------------------------------"
    1103             :             END IF
    1104             :          END IF
    1105             : 
    1106             :          CALL cp_print_key_finished_output(energy_unit, logger, &
    1107        1062 :                                            tddfpt_print_section, "DETAILED_ENERGY")
    1108             :          CALL cp_print_key_finished_output(log_unit, logger, &
    1109        1062 :                                            tddfpt_print_section, "ITERATION_INFO")
    1110             :       ELSE
    1111             :          CALL cp_warn(__LOCATION__, &
    1112           0 :                       "Skipping TDDFPT wavefunction optimization")
    1113             :       END IF
    1114             : 
    1115             :       IF (ASSOCIATED(matrix_ks_oep)) THEN
    1116             :          IF (tddfpt_control%dipole_form == tddfpt_dipole_velocity) THEN
    1117             :             CALL cp_warn(__LOCATION__, &
    1118             :                          "Transition dipole moments and oscillator strengths are likely to be incorrect "// &
    1119             :                          "when computed using an orbital energy correction XC-potential together with "// &
    1120             :                          "the velocity form of dipole transition integrals")
    1121             :          END IF
    1122             :       END IF
    1123             : 
    1124             :       ! *** print summary information ***
    1125        1062 :       log_unit = cp_logger_get_default_io_unit(logger)
    1126             : 
    1127             :       namd_print_section => section_vals_get_subs_vals( &
    1128             :                             tddfpt_print_section, &
    1129        1062 :                             "NAMD_PRINT")
    1130        1062 :       CALL section_vals_get(namd_print_section, explicit=explicit)
    1131        1062 :       IF (explicit) THEN
    1132             :          CALL tddfpt_write_newtonx_output(evects, &
    1133             :                                           evals, &
    1134             :                                           gs_mos, &
    1135             :                                           logger, &
    1136             :                                           tddfpt_print_section, &
    1137             :                                           matrix_s(1)%matrix, &
    1138             :                                           S_evects, &
    1139           2 :                                           sub_env)
    1140             :       END IF
    1141        3186 :       ALLOCATE (ostrength(nstates))
    1142        3936 :       ostrength = 0.0_dp
    1143             :       CALL tddfpt_print_summary(log_unit, &
    1144             :                                 evects, &
    1145             :                                 evals, &
    1146             :                                 ostrength, &
    1147             :                                 mult, &
    1148             :                                 dipole_op_mos_occ, &
    1149        1062 :                                 tddfpt_control%dipole_form)
    1150             :       CALL tddfpt_print_excitation_analysis( &
    1151             :          log_unit, &
    1152             :          evects, &
    1153             :          evals, &
    1154             :          gs_mos, &
    1155             :          matrix_s(1)%matrix, &
    1156        1062 :          min_amplitude=tddfpt_control%min_excitation_amplitude)
    1157             :       CALL tddfpt_print_nto_analysis(qs_env, &
    1158             :                                      evects, evals, &
    1159             :                                      ostrength, &
    1160             :                                      gs_mos, &
    1161             :                                      matrix_s(1)%matrix, &
    1162        1062 :                                      tddfpt_print_section)
    1163             : 
    1164        1062 :       IF (tddfpt_control%do_lrigpw) THEN
    1165             :          CALL lri_print_stat(qs_env, &
    1166             :                              ltddfpt=.TRUE., &
    1167          10 :                              tddfpt_lri_env=kernel_env%full_kernel%lri_env)
    1168             :       END IF
    1169             : 
    1170        1062 :       CALL timestop(handle)
    1171        4248 :    END SUBROUTINE tddfpt_energies
    1172             : 
    1173             : ! **************************************************************************************************
    1174             : !> \brief Perform singlet and triplet computations for subsequent TDDFPT-SOC calculation.
    1175             : !> \param qs_env  Quickstep environment
    1176             : !> \param nstates number of requested  exited states
    1177             : !> \param work_matrices ...
    1178             : !> \param tddfpt_control ...
    1179             : !> \param logger ...
    1180             : !> \param tddfpt_print_section ...
    1181             : !> \param evects Eigenvector of the requested multiplicity
    1182             : !> \param evals Eigenvalue of the requested multiplicity
    1183             : !> \param ostrength Oscillatorstrength
    1184             : !> \param gs_mos ...
    1185             : !> \param tddfpt_section ...
    1186             : !> \param S_evects ...
    1187             : !> \param matrix_s ...
    1188             : !> \param kernel_env ...
    1189             : !> \param matrix_ks ...
    1190             : !> \param sub_env ...
    1191             : !> \param dipole_op_mos_occ ...
    1192             : !> \param lmult_tmp ...
    1193             : !> \param xc_section ...
    1194             : !> \param full_kernel_env ...
    1195             : !> \param kernel_env_admm_aux ...
    1196             : !> \par History
    1197             : !>    * 02.2023 created [Jan-Robert Vogt]
    1198             : !> \note Based on tddfpt2_methods and xas_tdp_utils.
    1199             : !> \note only the values of one multiplicity will be passed back for force calculations!
    1200             : ! **************************************************************************************************
    1201             : 
    1202           8 :    SUBROUTINE tddfpt_soc_energies(qs_env, nstates, work_matrices, &
    1203             :                                   tddfpt_control, logger, tddfpt_print_section, &
    1204             :                                   evects, evals, ostrength, &
    1205             :                                   gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
    1206             :                                   sub_env, dipole_op_mos_occ, lmult_tmp, xc_section, full_kernel_env, &
    1207             :                                   kernel_env_admm_aux)
    1208             : 
    1209             :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
    1210             :       INTEGER, INTENT(in)                                :: nstates
    1211             :       TYPE(tddfpt_work_matrices)                         :: work_matrices
    1212             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
    1213             :       TYPE(cp_logger_type), POINTER                      :: logger
    1214             :       TYPE(section_vals_type), POINTER                   :: tddfpt_print_section
    1215             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: evects
    1216             :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: evals, ostrength
    1217             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1218             :          POINTER                                         :: gs_mos
    1219             :       TYPE(section_vals_type), POINTER                   :: tddfpt_section
    1220             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: S_evects
    1221             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
    1222             :       TYPE(kernel_env_type)                              :: kernel_env
    1223             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
    1224             :       TYPE(tddfpt_subgroup_env_type)                     :: sub_env
    1225             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: dipole_op_mos_occ
    1226             :       LOGICAL, INTENT(in)                                :: lmult_tmp
    1227             :       TYPE(section_vals_type), POINTER                   :: xc_section
    1228             :       TYPE(full_kernel_env_type), TARGET                 :: full_kernel_env, kernel_env_admm_aux
    1229             : 
    1230             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_soc_energies'
    1231             : 
    1232             :       INTEGER                                            :: handle, ispin, istate, log_unit, mult, &
    1233             :                                                             nspins
    1234           8 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: evals_mult, ostrength_mult
    1235           8 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: evects_mult
    1236             : 
    1237           8 :       CALL timeset(routineN, handle)
    1238             : 
    1239             :       log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
    1240             :                                       "PROGRAM_BANNER", &
    1241           8 :                                       extension=".tddfptLog")
    1242           8 :       CALL tddfpt_soc_header(log_unit)
    1243             : 
    1244           8 :       nspins = SIZE(gs_mos)
    1245          76 :       ALLOCATE (evects_mult(nspins, nstates))
    1246          24 :       ALLOCATE (evals_mult(nstates))
    1247             : 
    1248             :       ! First multiplicity
    1249           8 :       IF (lmult_tmp) THEN
    1250           2 :          IF (log_unit > 0) THEN
    1251           1 :             WRITE (log_unit, "(1X,A)") "", &
    1252           1 :                "-------------------------------------------------------------------------------", &
    1253           1 :                "-                      TDDFPT SINGLET ENERGIES                                 -", &
    1254           2 :                "-------------------------------------------------------------------------------"
    1255             :          END IF
    1256           2 :          mult = 1
    1257             :       ELSE
    1258           6 :          IF (log_unit > 0) THEN
    1259           3 :             WRITE (log_unit, "(1X,A)") "", &
    1260           3 :                "-------------------------------------------------------------------------------", &
    1261           3 :                "-                      TDDFPT TRIPLET ENERGIES                                 -", &
    1262           6 :                "-------------------------------------------------------------------------------"
    1263             :          END IF
    1264           6 :          mult = 3
    1265             :       END IF
    1266             : 
    1267             :       CALL tddfpt_energies(qs_env, nstates, work_matrices, tddfpt_control, logger, &
    1268             :                            tddfpt_print_section, evects_mult, evals_mult, &
    1269             :                            gs_mos, tddfpt_section, S_evects, matrix_s, &
    1270             :                            kernel_env, matrix_ks, sub_env, ostrength_mult, &
    1271             :                            dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
    1272           8 :                            kernel_env_admm_aux)
    1273             : 
    1274             :       ! Clean up in between for full kernel
    1275           8 :       IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
    1276           8 :          IF (tddfpt_control%do_admm) CALL release_kernel_env(kernel_env%admm_kernel)
    1277           8 :          CALL release_kernel_env(kernel_env%full_kernel)
    1278           8 :          CALL tddfpt_release_work_matrices(work_matrices, sub_env)
    1279             :          CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, tddfpt_control%do_hfx, &
    1280             :                                           tddfpt_control%do_admm, tddfpt_control%do_hfxlr, &
    1281           8 :                                           tddfpt_control%do_exck, qs_env, sub_env)
    1282             :       END IF
    1283             : 
    1284          30 :       DO istate = 1, nstates
    1285          52 :          DO ispin = 1, nspins
    1286          44 :             CALL cp_fm_release(S_evects(ispin, istate))
    1287             :          END DO
    1288             :       END DO
    1289             : 
    1290          30 :       DO istate = 1, nstates
    1291          52 :          DO ispin = 1, nspins
    1292             :             CALL fm_pool_create_fm( &
    1293             :                work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
    1294          44 :                S_evects(ispin, istate))
    1295             :          END DO
    1296             :       END DO
    1297             : 
    1298           8 :       tddfpt_control%rks_triplets = lmult_tmp
    1299             : 
    1300             :       ! Second multiplicity
    1301           8 :       IF (lmult_tmp) THEN
    1302           2 :          IF (log_unit > 0) THEN
    1303           1 :             WRITE (log_unit, "(1X,A)") "", &
    1304           1 :                "                      singlet excitations finished                             ", &
    1305           1 :                "                                                                               ", &
    1306           1 :                "-------------------------------------------------------------------------------", &
    1307           1 :                "-                      TDDFPT TRIPLET ENERGIES                                -", &
    1308           2 :                "-------------------------------------------------------------------------------"
    1309             :          END IF !log_unit
    1310           2 :          mult = 3
    1311             :       ELSE
    1312           6 :          IF (log_unit > 0) THEN
    1313           3 :             WRITE (log_unit, "(1X,A)") "", &
    1314           3 :                "                      triplet excitations finished                             ", &
    1315           3 :                "                                                                               ", &
    1316           3 :                "-------------------------------------------------------------------------------", &
    1317           3 :                "-                      TDDFPT SINGLET ENERGIES                                -", &
    1318           6 :                "-------------------------------------------------------------------------------"
    1319             :          END IF !log_unit
    1320           6 :          mult = 1
    1321             :       END IF
    1322             : 
    1323             :       CALL tddfpt_energies(qs_env, nstates, work_matrices, tddfpt_control, logger, &
    1324             :                            tddfpt_print_section, evects, evals, &
    1325             :                            gs_mos, tddfpt_section, S_evects, matrix_s, &
    1326             :                            kernel_env, matrix_ks, sub_env, ostrength, &
    1327             :                            dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
    1328           8 :                            kernel_env_admm_aux)
    1329             : 
    1330             :       ! Compute perturbative SOC correction
    1331             :       ! Order should always be singlet triplet in tddfpt_soc
    1332           8 :       IF (lmult_tmp) THEN
    1333           2 :          CALL tddfpt_soc(qs_env, evals_mult, evals, evects_mult, evects, gs_mos) !mult=singlet
    1334             :       ELSE
    1335           6 :          CALL tddfpt_soc(qs_env, evals, evals_mult, evects, evects_mult, gs_mos) !mult=triplet
    1336             :       END IF
    1337             : 
    1338             :       ! deallocate the additional multiplicity
    1339          16 :       DO ispin = 1, SIZE(evects_mult, 1)
    1340          38 :          DO istate = 1, SIZE(evects_mult, 2)
    1341          30 :             CALL cp_fm_release(evects_mult(ispin, istate))
    1342             :          END DO
    1343             :       END DO
    1344           8 :       DEALLOCATE (evects_mult, evals_mult, ostrength_mult)
    1345             : 
    1346           8 :       CALL timestop(handle)
    1347             : 
    1348          16 :    END SUBROUTINE tddfpt_soc_energies
    1349             : 
    1350             : END MODULE qs_tddfpt2_methods

Generated by: LCOV version 1.15