LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_soc.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:32ddf85) Lines: 606 620 97.7 %
Date: 2025-05-17 08:08:58 Functions: 6 7 85.7 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : MODULE qs_tddfpt2_soc
       8             :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
       9             :                                               gto_basis_set_type
      10             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      11             :    USE cp_cfm_diag,                     ONLY: cp_cfm_heevd
      12             :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      13             :                                               cp_cfm_get_info,&
      14             :                                               cp_cfm_get_submatrix,&
      15             :                                               cp_cfm_release,&
      16             :                                               cp_cfm_type,&
      17             :                                               cp_fm_to_cfm
      18             :    USE cp_control_types,                ONLY: dft_control_type,&
      19             :                                               qs_control_type,&
      20             :                                               tddfpt2_control_type
      21             :    USE cp_dbcsr_api,                    ONLY: &
      22             :         dbcsr_copy, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, &
      23             :         dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_get_info, dbcsr_multiply, &
      24             :         dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_type, dbcsr_type_no_symmetry
      25             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_print,&
      26             :                                               dbcsr_reserve_all_blocks
      27             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      28             :                                               copy_fm_to_dbcsr,&
      29             :                                               cp_dbcsr_sm_fm_multiply
      30             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale,&
      31             :                                               cp_fm_transpose,&
      32             :                                               cp_fm_uplo_to_full
      33             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      34             :                                               cp_fm_struct_release,&
      35             :                                               cp_fm_struct_type
      36             :    USE cp_fm_types,                     ONLY: &
      37             :         cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, &
      38             :         cp_fm_set_all, cp_fm_set_element, cp_fm_to_fm_submat, cp_fm_type
      39             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      40             :                                               cp_logger_get_default_unit_nr,&
      41             :                                               cp_logger_type
      42             :    USE cp_output_handling,              ONLY: cp_print_key_unit_nr
      43             :    USE input_constants,                 ONLY: tddfpt_dipole_length
      44             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      45             :                                               section_vals_type,&
      46             :                                               section_vals_val_get
      47             :    USE kinds,                           ONLY: default_string_length,&
      48             :                                               dp
      49             :    USE lebedev,                         ONLY: deallocate_lebedev_grids,&
      50             :                                               get_number_of_lebedev_grid,&
      51             :                                               init_lebedev_grids,&
      52             :                                               lebedev_grid
      53             :    USE mathlib,                         ONLY: get_diag
      54             :    USE memory_utilities,                ONLY: reallocate
      55             :    USE message_passing,                 ONLY: mp_para_env_type
      56             :    USE orbital_pointers,                ONLY: indso,&
      57             :                                               nsoset
      58             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      59             :    USE particle_types,                  ONLY: particle_type
      60             :    USE physcon,                         ONLY: evolt,&
      61             :                                               wavenumbers
      62             :    USE qs_environment_types,            ONLY: get_qs_env,&
      63             :                                               qs_environment_type
      64             :    USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
      65             :                                               create_grid_atom,&
      66             :                                               grid_atom_type
      67             :    USE qs_harmonics_atom,               ONLY: allocate_harmonics_atom,&
      68             :                                               create_harmonics_atom,&
      69             :                                               get_maxl_CG,&
      70             :                                               harmonics_atom_type
      71             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      72             :                                               get_qs_kind_set,&
      73             :                                               qs_kind_type
      74             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      75             :                                               mo_set_type
      76             :    USE qs_tddfpt2_soc_types,            ONLY: soc_atom_create,&
      77             :                                               soc_atom_env_type,&
      78             :                                               soc_atom_release,&
      79             :                                               soc_env_create,&
      80             :                                               soc_env_release,&
      81             :                                               soc_env_type
      82             :    USE qs_tddfpt2_soc_utils,            ONLY: dip_vel_op,&
      83             :                                               resort_evects,&
      84             :                                               soc_contract_evect,&
      85             :                                               soc_dipole_operator
      86             :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos
      87             :    USE soc_pseudopotential_methods,     ONLY: V_SOC_xyz_from_pseudopotential
      88             :    USE spherical_harmonics,             ONLY: clebsch_gordon,&
      89             :                                               clebsch_gordon_deallocate,&
      90             :                                               clebsch_gordon_init
      91             :    USE xas_tdp_atom,                    ONLY: compute_sphi_so,&
      92             :                                               integrate_soc_atoms,&
      93             :                                               truncate_radial_grid
      94             :    USE xas_tdp_utils,                   ONLY: rcs_amew_soc_elements
      95             : #include "./base/base_uses.f90"
      96             : 
      97             :    IMPLICIT NONE
      98             : 
      99             :    PRIVATE
     100             : 
     101             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_soc'
     102             : 
     103             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
     104             : 
     105             :    INTEGER, PARAMETER, PRIVATE          :: nderivs = 3
     106             :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
     107             : 
     108             :    !A helper type for SOC
     109             :    TYPE dbcsr_soc_package_type
     110             :       TYPE(dbcsr_type), POINTER     :: dbcsr_sg => Null()
     111             :       TYPE(dbcsr_type), POINTER     :: dbcsr_tp => Null()
     112             :       TYPE(dbcsr_type), POINTER     :: dbcsr_sc => Null()
     113             :       TYPE(dbcsr_type), POINTER     :: dbcsr_sf => Null()
     114             :       TYPE(dbcsr_type), POINTER     :: dbcsr_prod => Null()
     115             :       TYPE(dbcsr_type), POINTER     :: dbcsr_ovlp => Null()
     116             :       TYPE(dbcsr_type), POINTER     :: dbcsr_tmp => Null()
     117             :       TYPE(dbcsr_type), POINTER     :: dbcsr_work => Null()
     118             :    END TYPE dbcsr_soc_package_type
     119             : 
     120             :    PUBLIC :: tddfpt_soc
     121             : 
     122             : ! **************************************************************************************************
     123             : 
     124             : CONTAINS
     125             : 
     126             : ! **************************************************************************************************
     127             : !> \brief Perform TDDFPT-SOC calculation.
     128             : !> \param qs_env  Quickstep environment
     129             : !> \param evals_a eigenvalues for the singlet states
     130             : !> \param evals_b eigenvalues for the triplet states
     131             : !> \param evects_a eigenvectors for the singlet states
     132             : !> \param evects_b eigenvectors for the triplet states
     133             : !> \param gs_mos ground state orbitlas from the TDDFPT calculation
     134             : !> \par History
     135             : !>    * 02.2023 created [Jan-Robert Vogt]
     136             : !> \note Based on tddfpt2_methods and xas_tdp_utils.
     137             : !> \note Only rcs for now, but written with the addition of os in mind
     138             : ! **************************************************************************************************
     139             : 
     140           8 :    SUBROUTINE tddfpt_soc(qs_env, evals_a, evals_b, evects_a, evects_b, gs_mos)
     141             : 
     142             :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
     143             :       REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: evals_a, evals_b
     144             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects_a, evects_b
     145             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     146             :          POINTER                                         :: gs_mos
     147             : 
     148             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_soc'
     149             : 
     150             :       INTEGER                                            :: handle, istate, log_unit
     151             :       LOGICAL                                            :: do_os
     152             :       TYPE(cp_logger_type), POINTER                      :: logger
     153             :       TYPE(dft_control_type), POINTER                    :: dft_control
     154             : 
     155           8 :       CALL timeset(routineN, handle)
     156             : 
     157           8 :       logger => cp_get_default_logger()
     158           8 :       IF (logger%para_env%is_source()) THEN
     159           4 :          log_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     160             :       ELSE
     161           4 :          log_unit = -1
     162             :       END IF
     163             : 
     164           8 :       IF (log_unit > 0) THEN
     165           4 :          WRITE (log_unit, "(1X,A)") "", &
     166           4 :             "-------------------------------------------------------------------------------", &
     167           4 :             "-                         START SOC CALCULATIONS                              -", &
     168           8 :             "-------------------------------------------------------------------------------"
     169             :       END IF
     170             : 
     171           8 :       NULLIFY (dft_control)
     172           8 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     173           8 :       do_os = dft_control%uks .OR. dft_control%roks
     174           0 :       IF (do_os) CPABORT("SOC only implemented for closed shell.")
     175             : 
     176           8 :       IF (log_unit > 0) THEN
     177           4 :          WRITE (log_unit, '(A)') "Starting from TDDFPT Excited States:"
     178           4 :          WRITE (log_unit, '(A)') "      STATE          SINGLET/eV         TRIPLET/eV"
     179          15 :          DO istate = 1, SIZE(evals_a)
     180          15 :             WRITE (log_unit, '(6X,I3,11X,F10.5,6X,F10.5)') istate, evals_a(istate)*evolt, evals_b(istate)*evolt
     181             :          END DO
     182             :       END IF
     183             : 
     184           8 :       IF (log_unit > 0) WRITE (log_unit, '(A)') "Starting restricted closed shell:"
     185           8 :       CALL tddfpt_soc_rcs(qs_env, evals_a, evals_b, evects_a, evects_b, log_unit, gs_mos)
     186             : 
     187           8 :       IF (log_unit > 0) THEN
     188           4 :          WRITE (log_unit, '(A,/,A)') "SOC Calculation terminated", &
     189           8 :             "Returning to TDDFPT for Force calculation and deallocations"
     190             :       END IF
     191             : 
     192           8 :       CALL timestop(handle)
     193             : 
     194           8 :    END SUBROUTINE
     195             : 
     196             : ! **************************************************************************************************
     197             : !> \brief Will perform the soc-calculation for restricted-closed-shell systems
     198             : !> \param qs_env Quickstep Enviroment
     199             : !> \param evals_sing eigenvalues singlet states
     200             : !> \param evals_trip eigenvalues triplet states
     201             : !> \param evects_sing eigenvector singlet states
     202             : !> \param evects_trip eigenvectors triplet states
     203             : !> \param log_unit default log_unit for convinients
     204             : !> \param gs_mos ground state MOs from TDDFPT
     205             : !> \par History
     206             : !>      * created 02.2023 [Jan-Robert Vogt]
     207             : !> \note Mostly copied and modified from xas_tdp_utils include_rcs_soc
     208             : ! **************************************************************************************************
     209           8 :    SUBROUTINE tddfpt_soc_rcs(qs_env, evals_sing, evals_trip, evects_sing, evects_trip, log_unit, gs_mos)
     210             : 
     211             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     212             :       REAL(kind=dp), DIMENSION(:), INTENT(IN), TARGET    :: evals_sing, evals_trip
     213             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects_sing, evects_trip
     214             :       INTEGER                                            :: log_unit
     215             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     216             :          POINTER                                         :: gs_mos
     217             : 
     218             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tddfpt_soc_rcs'
     219             : 
     220             :       CHARACTER(len=3)                                   :: mult
     221             :       INTEGER                                            :: dipole_form, group, handle, iex, ii, &
     222             :                                                             isg, istate, itp, jj, nactive, nao, &
     223             :                                                             nex, npcols, nprows, nsg, ntot, ntp
     224           8 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: evects_sort
     225           8 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, col_dist, row_blk_size, &
     226           8 :                                                             row_dist, row_dist_new
     227           8 :       INTEGER, DIMENSION(:, :), POINTER                  :: pgrid
     228             :       LOGICAL                                            :: print_ev, print_some, print_splitting, &
     229             :                                                             print_wn
     230             :       REAL(dp)                                           :: eps_filter, soc_gst, sqrt2, unit_scale
     231           8 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: diag, tmp_evals
     232           8 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: img_tmp, real_tmp
     233           8 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: gstp_block, mo_soc_x, mo_soc_y, mo_soc_z
     234             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     235             :       TYPE(cp_cfm_type)                                  :: evecs_cfm, hami_cfm
     236             :       TYPE(cp_fm_struct_type), POINTER                   :: full_struct, gstp_struct, prod_struct, &
     237             :                                                             vec_struct, work_struct
     238             :       TYPE(cp_fm_type)                                   :: gstp_fm, img_fm, prod_fm, real_fm, &
     239             :                                                             tmp_fm, vec_soc_x, vec_soc_y, &
     240             :                                                             vec_soc_z, work_fm
     241             :       TYPE(cp_fm_type), POINTER                          :: gs_coeffs, tp_coeffs
     242             :       TYPE(cp_logger_type), POINTER                      :: logger
     243             :       TYPE(dbcsr_distribution_type), POINTER             :: coeffs_dist, dbcsr_dist, prod_dist
     244           8 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     245             :       TYPE(dbcsr_soc_package_type)                       :: dbcsr_soc_package
     246             :       TYPE(dbcsr_type), POINTER                          :: dbcsr_dummy, dbcsr_ovlp, dbcsr_prod, &
     247             :                                                             dbcsr_sg, dbcsr_tmp, dbcsr_tp, &
     248             :                                                             dbcsr_work, orb_soc_x, orb_soc_y, &
     249             :                                                             orb_soc_z
     250           8 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     251             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     252             :       TYPE(section_vals_type), POINTER                   :: soc_print_section, soc_section, &
     253             :                                                             tddfpt_print_section
     254             :       TYPE(soc_atom_env_type), POINTER                   :: soc_atom_env
     255           8 :       TYPE(soc_env_type), TARGET                         :: soc_env
     256             : 
     257           8 :       CALL timeset(routineN, handle)
     258             : 
     259           8 :       NULLIFY (logger, pgrid, row_dist, row_dist_new)
     260           8 :       NULLIFY (col_dist, col_blk_size, row_blk_size, matrix_s)
     261           8 :       NULLIFY (gstp_struct, tddfpt_print_section, soc_print_section, soc_section)
     262             : 
     263           8 :       logger => cp_get_default_logger()
     264           8 :       tddfpt_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
     265           8 :       soc_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT%SOC_PRINT")
     266           8 :       soc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%SOC")
     267             : 
     268           8 :       CALL section_vals_val_get(soc_section, "EPS_FILTER", r_val=eps_filter)
     269             : 
     270           8 :       nsg = SIZE(evals_sing)   ! Number of excited singlet states
     271           8 :       ntp = SIZE(evals_trip)   ! Number of excited triplet states
     272           8 :       nex = nsg                ! Number of excited states of each multiplicity
     273           8 :       ntot = 1 + nsg + 3*ntp ! Number of (GS + S + T^-1 + T^0 + T^1)
     274             : 
     275             :       ! Initzialize Working environment
     276             :       CALL inititialize_soc(qs_env, soc_atom_env, soc_env, &
     277           8 :                             evects_sing, evects_trip, dipole_form, gs_mos)
     278             : 
     279           8 :       NULLIFY (mos)
     280           8 :       CALL get_qs_env(qs_env, mos=mos, para_env=para_env, blacs_env=blacs_env)
     281           8 :       CALL get_mo_set(mos(1), nao=nao, homo=nactive)
     282             : 
     283             :       ! this will create the H^SOC in an atomic basis
     284           8 :       CALL integrate_soc_atoms(soc_env%orb_soc, qs_env=qs_env, soc_atom_env=soc_atom_env)
     285           8 :       CALL soc_atom_release(soc_atom_env)
     286             : 
     287             :       !! Point at H^SOC and MOs for better readablity
     288           8 :       NULLIFY (tp_coeffs, orb_soc_x, orb_soc_y, orb_soc_z)
     289           8 :       tp_coeffs => soc_env%b_coeff
     290           8 :       soc_env%evals_a => evals_sing
     291           8 :       soc_env%evals_b => evals_trip
     292           8 :       orb_soc_x => soc_env%orb_soc(1)%matrix
     293           8 :       orb_soc_y => soc_env%orb_soc(2)%matrix
     294           8 :       orb_soc_z => soc_env%orb_soc(3)%matrix
     295             : 
     296             :       !! Create a matrix-structure, which links all states in this calculation
     297           8 :       NULLIFY (full_struct)
     298           8 :       CALL cp_fm_struct_create(full_struct, context=blacs_env, para_env=para_env, nrow_global=ntot, ncol_global=ntot)
     299           8 :       CALL cp_fm_create(real_fm, full_struct)
     300           8 :       CALL cp_fm_create(img_fm, full_struct)
     301           8 :       CALL cp_fm_set_all(real_fm, 0.0_dp)
     302           8 :       CALL cp_fm_set_all(img_fm, 0.0_dp)
     303             : 
     304             :       !  Put the excitation energies on the diagonal of the real matrix
     305          30 :       DO isg = 1, nsg
     306          30 :          CALL cp_fm_set_element(real_fm, 1 + isg, 1 + isg, evals_sing(isg))
     307             :       END DO
     308          30 :       DO itp = 1, ntp
     309             :          ! first T^-1, then T^0, then T^+1
     310          22 :          CALL cp_fm_set_element(real_fm, 1 + itp + nsg, 1 + itp + nsg, evals_trip(itp))
     311          22 :          CALL cp_fm_set_element(real_fm, 1 + itp + ntp + nsg, 1 + itp + ntp + nsg, evals_trip(itp))
     312          30 :          CALL cp_fm_set_element(real_fm, 1 + itp + 2*ntp + nsg, 1 + itp + 2*ntp + nsg, evals_trip(itp))
     313             :       END DO
     314             : 
     315             :       !!Create the dbcsr structures for this calculations
     316           8 :       CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist)
     317             :       CALL dbcsr_distribution_get(dbcsr_dist, group=group, row_dist=row_dist, pgrid=pgrid, &
     318           8 :                                   npcols=npcols, nprows=nprows)
     319          32 :       ALLOCATE (col_dist(nex), row_dist_new(nex))                   ! Split for each excitation
     320          30 :       DO iex = 1, nex
     321          22 :          col_dist(iex) = MODULO(npcols - iex, npcols)
     322          30 :          row_dist_new(iex) = MODULO(nprows - iex, nprows)
     323             :       END DO
     324           8 :       ALLOCATE (coeffs_dist, prod_dist)
     325             :       CALL dbcsr_distribution_new(coeffs_dist, group=group, pgrid=pgrid, row_dist=row_dist, &
     326           8 :                                   col_dist=col_dist)
     327             :       CALL dbcsr_distribution_new(prod_dist, group=group, pgrid=pgrid, row_dist=row_dist_new, &
     328           8 :                                   col_dist=col_dist)
     329             : 
     330             :       !! Create the matrices
     331          16 :       ALLOCATE (col_blk_size(nex))
     332          30 :       col_blk_size = nactive
     333             : 
     334           8 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     335           8 :       CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=row_blk_size)
     336             : 
     337             :       !! The Eigenvectors for Sg und Tp will be dived into their diffrent components again
     338           8 :       ALLOCATE (dbcsr_sg, dbcsr_tp, dbcsr_work, dbcsr_ovlp, dbcsr_tmp, dbcsr_prod)
     339             :       CALL dbcsr_create(matrix=dbcsr_sg, name="SINGLETS", matrix_type=dbcsr_type_no_symmetry, &
     340           8 :                         dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
     341             :       CALL dbcsr_create(matrix=dbcsr_tp, name="TRIPLETS", matrix_type=dbcsr_type_no_symmetry, &
     342           8 :                         dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
     343             :       CALL dbcsr_create(matrix=dbcsr_work, name="WORK", matrix_type=dbcsr_type_no_symmetry, &
     344           8 :                         dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
     345             :       CALL dbcsr_create(matrix=dbcsr_prod, name="PROD", matrix_type=dbcsr_type_no_symmetry, &
     346           8 :                         dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
     347             :       CALL dbcsr_create(matrix=dbcsr_ovlp, name="OVLP", matrix_type=dbcsr_type_no_symmetry, &
     348           8 :                         dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
     349             : 
     350          30 :       col_blk_size = 1
     351             :       CALL dbcsr_create(matrix=dbcsr_tmp, name="TMP", matrix_type=dbcsr_type_no_symmetry, &
     352           8 :                         dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
     353           8 :       CALL dbcsr_reserve_all_blocks(dbcsr_tmp)
     354             : 
     355             :       IF (debug_this_module) THEN
     356             :          ALLOCATE (dbcsr_dummy)
     357             :          CALL dbcsr_create(matrix=dbcsr_dummy, name="DUMMY", matrix_type=dbcsr_type_no_symmetry, &
     358             :                            dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
     359             :          CALL dbcsr_reserve_all_blocks(dbcsr_dummy)
     360             :       END IF
     361             : 
     362             :       !! This work dbcsr matrix will be packed together for easy transfer to other subroutines
     363           8 :       dbcsr_soc_package%dbcsr_sg => dbcsr_sg
     364           8 :       dbcsr_soc_package%dbcsr_tp => dbcsr_tp
     365           8 :       dbcsr_soc_package%dbcsr_work => dbcsr_work
     366           8 :       dbcsr_soc_package%dbcsr_ovlp => dbcsr_ovlp
     367           8 :       dbcsr_soc_package%dbcsr_prod => dbcsr_prod
     368           8 :       dbcsr_soc_package%dbcsr_tmp => dbcsr_tmp
     369             : 
     370             :       !Filling the coeffs matrices by copying from the stored fms
     371           8 :       CALL copy_fm_to_dbcsr(soc_env%a_coeff, dbcsr_sg)
     372           8 :       CALL copy_fm_to_dbcsr(soc_env%b_coeff, dbcsr_tp)
     373             : 
     374             :       !Create the work and helper fms
     375             :       !CALL get_mo_set(mos(1),mo_coeff=gs_coeffs)
     376           8 :       gs_coeffs => gs_mos(1)%mos_occ
     377           8 :       CALL cp_fm_get_info(gs_coeffs, matrix_struct=vec_struct)
     378           8 :       CALL cp_fm_create(vec_soc_x, vec_struct)
     379           8 :       CALL cp_fm_create(vec_soc_y, vec_struct)
     380           8 :       CALL cp_fm_create(vec_soc_z, vec_struct)
     381             : 
     382             :       CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
     383           8 :                                nrow_global=nactive, ncol_global=nactive)
     384           8 :       CALL cp_fm_create(prod_fm, prod_struct)
     385             : 
     386             :       CALL cp_fm_struct_create(work_struct, context=blacs_env, para_env=para_env, &
     387           8 :                                nrow_global=nex, ncol_global=nex)
     388           8 :       CALL cp_fm_create(work_fm, work_struct)
     389           8 :       CALL cp_fm_create(tmp_fm, work_struct)
     390             : 
     391             :       IF (log_unit > 0 .AND. debug_this_module) THEN
     392             :          WRITE (log_unit, '(A)') "Starting Precomputations"
     393             :       END IF
     394             : 
     395             :       !! Begin with the precomputation
     396             :       !! Prefactor due to rcs.
     397             :       !! The excitation is presented as a linear combination of
     398             :       !! alpha and beta, where dalpha=-dbeta for triplet exc.
     399             :       !! and dalpha=dbeta for the singlet case
     400           8 :       sqrt2 = SQRT(2.0_dp)
     401             : 
     402             :       !! Precompute the <phi_i^0|H^SOC|phi_j^0> matrix elements
     403          24 :       ALLOCATE (diag(nactive))
     404          64 :       ALLOCATE (mo_soc_x(nactive, nactive), mo_soc_y(nactive, nactive), mo_soc_z(nactive, nactive))
     405             : 
     406             :       !! This will be the GS|H|GS contribution needed for all couplings
     407           8 :       CALL cp_dbcsr_sm_fm_multiply(orb_soc_x, gs_coeffs, vec_soc_x, ncol=nactive)
     408             :       CALL parallel_gemm('T', 'N', nactive, &
     409             :                          nactive, nao, 1.0_dp, &
     410             :                          gs_coeffs, vec_soc_x, &
     411           8 :                          0.0_dp, prod_fm)
     412           8 :       CALL cp_fm_get_submatrix(prod_fm, mo_soc_x)
     413             : 
     414           8 :       CALL cp_dbcsr_sm_fm_multiply(orb_soc_y, gs_coeffs, vec_soc_y, ncol=nactive)
     415             :       CALL parallel_gemm('T', 'N', nactive, &
     416             :                          nactive, nao, 1.0_dp, &
     417             :                          gs_coeffs, vec_soc_y, &
     418           8 :                          0.0_dp, prod_fm)
     419           8 :       CALL cp_fm_get_submatrix(prod_fm, mo_soc_y)
     420             : 
     421           8 :       CALL cp_dbcsr_sm_fm_multiply(orb_soc_z, gs_coeffs, vec_soc_z, ncol=nactive)
     422             :       CALL parallel_gemm('T', 'N', nactive, &
     423             :                          nactive, nao, 1.0_dp, &
     424             :                          gs_coeffs, vec_soc_z, &
     425           8 :                          0.0_dp, prod_fm)
     426           8 :       CALL cp_fm_get_submatrix(prod_fm, mo_soc_z)
     427             : 
     428             :       !  Only have SOC between singlet-triplet triplet-triplet and ground_state-triplet, the resulting
     429             :       !  matrix is Hermitian i.e. the real part is symmetric and the imaginary part is anti-symmetric.
     430             :       !  Can only fill upper half
     431             :       IF (log_unit > 0 .AND. debug_this_module) THEN
     432             :          WRITE (log_unit, '(A)') "Starting Ground-State contributions..."
     433             :       END IF
     434             :       !Start with the ground state/triplet SOC, SOC*gs_coeffs already computed above
     435             :       !note: we are computing <0|H|T>, but have SOC*gs_coeffs instead of gs_coeffs*SOC in store.
     436             :       !Since the SOC Hamiltonian is anti-symmetric, a - signs pops up in the gemms below
     437             :       CALL cp_fm_struct_create(gstp_struct, context=blacs_env, para_env=para_env, &
     438           8 :                                nrow_global=ntp*nactive, ncol_global=nactive)
     439           8 :       CALL cp_fm_create(gstp_fm, gstp_struct)
     440          24 :       ALLOCATE (gstp_block(nactive, nactive))
     441             : 
     442             :       !gs-triplet with Ms=+-1, imaginary part
     443             :       ! <T+-1|H_x|GS>
     444             :       ! -1 to change to <GS|H|T+-1>
     445             :       CALL parallel_gemm('T', 'N', nactive*ntp, &
     446             :                          nactive, nao, -1.0_dp, &
     447             :                          tp_coeffs, vec_soc_x, &
     448           8 :                          0.0_dp, gstp_fm)
     449             : 
     450             :       !! Seperate them into the different states again (nactive x nactive)
     451          30 :       DO itp = 1, ntp
     452             :          CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, &
     453             :                                   start_row=(itp - 1)*nactive + 1, start_col=1, &
     454          22 :                                   n_rows=nactive, n_cols=nactive)
     455          22 :          diag(:) = get_diag(gstp_block)
     456         176 :          soc_gst = SUM(diag)
     457             :          !!  <0|H_x|T^-1>
     458          22 :          CALL cp_fm_set_element(img_fm, 1, 1 + nsg + itp, -1_dp*soc_gst)
     459             :          !! <0|H_x|T^+1>
     460          30 :          CALL cp_fm_set_element(img_fm, 1, 1 + nsg + 2*ntp + itp, soc_gst)
     461             :       END DO
     462             : 
     463             :       IF (debug_this_module .AND. (log_unit > 0)) THEN
     464             :          WRITE (log_unit, "(A,2F18.6)") "<0|H_x|T^-1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
     465             :          WRITE (log_unit, "(A,2F18.6)") "<0|H_x|T^+1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
     466             :       END IF
     467             : 
     468             :       !gs-triplet with Ms=+-1, real part
     469             :       ! <T+-1|H_y|GS>
     470             :       CALL parallel_gemm('T', 'N', nactive*ntp, &
     471             :                          nactive, nao, -1.0_dp, &
     472             :                          tp_coeffs, vec_soc_y, &
     473           8 :                          0.0_dp, gstp_fm)
     474             : 
     475          30 :       DO itp = 1, ntp
     476             :          CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, start_row=(itp - 1)*nactive + 1, &
     477          22 :                                   start_col=1, n_rows=nactive, n_cols=nactive)
     478          22 :          diag(:) = get_diag(gstp_block)
     479         176 :          soc_gst = SUM(diag)
     480             :          ! <0|H_y|T^-1>
     481          22 :          CALL cp_fm_set_element(real_fm, 1, 1 + nsg + itp, -1.0_dp*soc_gst)
     482             :          ! <0|H_y|T^+1>
     483          30 :          CALL cp_fm_set_element(real_fm, 1, 1 + nsg + 2*ntp + itp, -1.0_dp*soc_gst)
     484             :       END DO
     485             : 
     486             :       IF (debug_this_module .AND. (log_unit > 0)) THEN
     487             :          WRITE (log_unit, "(A,2F18.6)") "<0|H_y|T^-1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
     488             :          WRITE (log_unit, "(A,2F18.6)") "<0|H_y|T^+1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
     489             :       END IF
     490             : 
     491             :       !gs-triplet with Ms=0, purely imaginary
     492             :       !< T0|H_z|GS>
     493             :       CALL parallel_gemm('T', 'N', nactive*ntp, &
     494             :                          nactive, nao, -1.0_dp, &
     495             :                          tp_coeffs, vec_soc_z, &
     496           8 :                          0.0_dp, gstp_fm)
     497             : 
     498          30 :       DO itp = 1, ntp
     499             :          CALL cp_fm_get_submatrix(fm=gstp_fm, &
     500             :                                   target_m=gstp_block, &
     501             :                                   start_row=(itp - 1)*nactive + 1, &
     502             :                                   start_col=1, &
     503             :                                   n_rows=nactive, &
     504          22 :                                   n_cols=nactive)
     505          22 :          diag(:) = get_diag(gstp_block)
     506         176 :          soc_gst = sqrt2*SUM(diag)
     507          30 :          CALL cp_fm_set_element(img_fm, 1, 1 + nsg + ntp + itp, soc_gst)
     508             :       END DO
     509             : 
     510             :       IF (debug_this_module .AND. (log_unit > 0)) THEN
     511             :          WRITE (log_unit, "(A,2F18.6)") "<0|H_z|T> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
     512             :       END IF
     513             : 
     514             :       !! After all::
     515             :       !! T-1 :: -<0|H_x|T^-1> - <0|H_y|T^-1> at 1+nsg+itp
     516             :       !! T+1 :: <0|H_x|T^-1> - <0|H_y|T^-1> at 1+nsp+2tnp+itp
     517             :       !! T0  :: < T0|H_z|GS>                at 1+nsg+ntp+itp
     518             : 
     519             :       !gs clean-up
     520           8 :       CALL cp_fm_release(prod_fm)
     521           8 :       CALL cp_fm_release(vec_soc_x)
     522           8 :       CALL cp_fm_release(vec_soc_y)
     523           8 :       CALL cp_fm_release(vec_soc_z)
     524           8 :       CALL cp_fm_release(gstp_fm)
     525           8 :       CALL cp_fm_struct_release(gstp_struct)
     526           8 :       CALL cp_fm_struct_release(prod_struct)
     527           8 :       DEALLOCATE (gstp_block)
     528             : 
     529             :       IF (log_unit > 0 .AND. debug_this_module) THEN
     530             :          WRITE (log_unit, '(A)') "Starting Singlet-Triplet contributions..."
     531             :       END IF
     532             : 
     533             :       !Now do the singlet-triplet SOC
     534             :       !start by computing the singlet-triplet overlap
     535             :       ! <S|T>
     536             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     537             :                           matrix_s(1)%matrix, &
     538             :                           dbcsr_tp, 0.0_dp, &
     539             :                           dbcsr_work, &
     540           8 :                           filter_eps=eps_filter)
     541             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     542             :                           dbcsr_sg, &
     543             :                           dbcsr_work, &
     544             :                           0.0_dp, dbcsr_ovlp, &
     545           8 :                           filter_eps=eps_filter)
     546             : 
     547             :       !singlet-triplet with Ms=+-1, imaginary part
     548             :       ! First precalculate <S|H_x|T+-1>
     549             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     550             :                           orb_soc_x, &
     551             :                           dbcsr_tp, 0.0_dp, &
     552             :                           dbcsr_work, &
     553           8 :                           filter_eps=eps_filter)
     554             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     555             :                           dbcsr_sg, &
     556             :                           dbcsr_work, 0.0_dp, &
     557             :                           dbcsr_prod, &
     558           8 :                           filter_eps=eps_filter)
     559             : 
     560             :       !! This will lead to:
     561             :       !! -1/sqrt(2)(<S|H_x|T> - <S|T> <GS|H_x|GS>)
     562             :       CALL rcs_amew_soc_elements(dbcsr_tmp, &
     563             :                                  dbcsr_prod, &
     564             :                                  dbcsr_ovlp, &
     565             :                                  mo_soc_x, &
     566             :                                  pref_trace=-1.0_dp, &
     567           8 :                                  pref_overall=-0.5_dp*sqrt2)
     568             : 
     569             :       !<S|H_x|T^-1>
     570             :       !! Convert to fm for transfer to img_fm
     571           8 :       CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
     572             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     573             :                               mtarget=img_fm, &
     574             :                               nrow=nex, &
     575             :                               ncol=nex, &
     576             :                               s_firstrow=1, &
     577             :                               s_firstcol=1, &
     578             :                               t_firstrow=2, &
     579           8 :                               t_firstcol=1 + nsg + 1)
     580             : 
     581             :       IF (debug_this_module) THEN
     582             :          WRITE (log_unit, "(/,A)") "<S|H_x|T^-1> in Hartree / cm^-1"
     583             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     584             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     585             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     586             :       END IF
     587             : 
     588             :       !<S|H_x|T^+1> takes a minus sign
     589           8 :       CALL cp_fm_scale(-1.0_dp, tmp_fm)
     590             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     591             :                               mtarget=img_fm, &
     592             :                               nrow=nex, &
     593             :                               ncol=nex, &
     594             :                               s_firstrow=1, &
     595             :                               s_firstcol=1, &
     596             :                               t_firstrow=2, &
     597           8 :                               t_firstcol=1 + nsg + 2*ntp + 1)
     598             : 
     599             :       IF (debug_this_module) THEN
     600             :          WRITE (log_unit, "(/,A)") "<S|H_x|T^+1> in Hartree / cm^-1"
     601             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     602             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     603             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     604             :       END IF
     605             : 
     606             :        !!singlet-triplet with Ms=+-1, real part
     607             :        !! Precompute <S|H_y|T>
     608             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     609             :                           orb_soc_y, dbcsr_tp, &
     610             :                           0.0_dp, dbcsr_work, &
     611           8 :                           filter_eps=eps_filter)
     612             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     613             :                           dbcsr_sg, dbcsr_work, &
     614             :                           0.0_dp, dbcsr_prod, &
     615           8 :                           filter_eps=eps_filter)
     616             : 
     617             :       !! This will lead to -1/sqrt(2)(<S|H_y|T> - <S|T> <GS|H_y|GS>)
     618             :       CALL rcs_amew_soc_elements(dbcsr_tmp, &
     619             :                                  dbcsr_prod, &
     620             :                                  dbcsr_ovlp, &
     621             :                                  mo_soc_y, &
     622             :                                  pref_trace=-1.0_dp, &
     623           8 :                                  pref_overall=-0.5_dp*sqrt2)
     624             : 
     625             :       !<S|H_y|T^-1>
     626           8 :       CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
     627             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     628             :                               mtarget=real_fm, &
     629             :                               nrow=nex, &
     630             :                               ncol=nex, &
     631             :                               s_firstrow=1, &
     632             :                               s_firstcol=1, &
     633             :                               t_firstrow=2, &
     634           8 :                               t_firstcol=1 + nsg + 1)
     635             : 
     636             :       IF (debug_this_module) THEN
     637             :          WRITE (log_unit, "(/,A)") "<S|H_y|T^-1> in Hartree / cm^-1"
     638             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     639             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     640             :          CALL dbcsr_print(dbcsr_tmp, unit_nr=log_unit)
     641             :       END IF
     642             : 
     643             :       !<S|H_y|T^+1>
     644             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     645             :                               mtarget=real_fm, &
     646             :                               nrow=nex, &
     647             :                               ncol=nex, &
     648             :                               s_firstrow=1, &
     649             :                               s_firstcol=1, &
     650             :                               t_firstrow=2, &
     651           8 :                               t_firstcol=1 + nsg + 2*ntp + 1)
     652             : 
     653             :       IF (debug_this_module) THEN
     654             :          WRITE (log_unit, "(/,A)") "<S|H_y|T^+1> in Hartree / cm^-1"
     655             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     656             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     657             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     658             :       END IF
     659             : 
     660             :       !singlet-triplet with Ms=0, purely imaginary
     661             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     662             :                           orb_soc_z, dbcsr_tp, &
     663             :                           0.0_dp, dbcsr_work, &
     664           8 :                           filter_eps=eps_filter)
     665             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     666             :                           dbcsr_sg, dbcsr_work, &
     667             :                           0.0_dp, dbcsr_prod, &
     668           8 :                           filter_eps=eps_filter)
     669             : 
     670             :       CALL rcs_amew_soc_elements(dbcsr_tmp, &
     671             :                                  dbcsr_prod, &
     672             :                                  dbcsr_ovlp, &
     673             :                                  mo_soc_z, &
     674             :                                  pref_trace=-1.0_dp, &
     675           8 :                                  pref_overall=1.0_dp)
     676             : 
     677             :       !<S|H_z|T^0>
     678           8 :       CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
     679             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     680             :                               mtarget=img_fm, &
     681             :                               nrow=nex, &
     682             :                               ncol=nex, &
     683             :                               s_firstrow=1, &
     684             :                               s_firstcol=1, &
     685             :                               t_firstrow=2, &
     686           8 :                               t_firstcol=1 + nsg + ntp + 1)
     687             : 
     688             :       IF (debug_this_module) THEN
     689             :          WRITE (log_unit, "(/,A)") "<S|H_z|T^0> in Hartree / cm^-1"
     690             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     691             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     692             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     693             :       END IF
     694             : 
     695             :       IF (log_unit > 0 .AND. debug_this_module) THEN
     696             :          WRITE (log_unit, '(A)') "Starting Triplet-Triplet contributions..."
     697             :       END IF
     698             : 
     699             :       !Now the triplet-triplet SOC
     700             :       !start by computing the overlap
     701             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     702             :                           matrix_s(1)%matrix, &
     703             :                           dbcsr_tp, 0.0_dp, &
     704             :                           dbcsr_work, &
     705           8 :                           filter_eps=eps_filter)
     706             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     707             :                           dbcsr_tp, dbcsr_work, &
     708             :                           0.0_dp, dbcsr_ovlp, &
     709           8 :                           filter_eps=eps_filter)
     710             : 
     711             :       !Ms=0 to Ms=+-1 SOC, imaginary part
     712             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     713             :                           orb_soc_x, dbcsr_tp, &
     714             :                           0.0_dp, dbcsr_work, &
     715           8 :                           filter_eps=eps_filter)
     716             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     717             :                           dbcsr_tp, dbcsr_work, &
     718             :                           0.0_dp, dbcsr_prod, &
     719           8 :                           filter_eps=eps_filter)
     720             : 
     721             :       CALL rcs_amew_soc_elements(dbcsr_tmp, &
     722             :                                  dbcsr_prod, &
     723             :                                  dbcsr_ovlp, &
     724             :                                  mo_soc_x, &
     725             :                                  pref_trace=1.0_dp, &
     726           8 :                                  pref_overall=-0.5_dp*sqrt2)
     727             : 
     728             :       !<T^0|H_x|T^+1>
     729           8 :       CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
     730             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     731             :                               mtarget=img_fm, &
     732             :                               nrow=nex, &
     733             :                               ncol=nex, &
     734             :                               s_firstrow=1, &
     735             :                               s_firstcol=1, &
     736             :                               t_firstrow=1 + nsg + ntp + 1, &
     737           8 :                               t_firstcol=1 + nsg + 2*ntp + 1)
     738             : 
     739             :       IF (debug_this_module) THEN
     740             :          WRITE (log_unit, "(/,A)") "<T^0|H_x|T^+1> in Hartree / cm^-1"
     741             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     742             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     743             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     744             :       END IF
     745             : 
     746             :       !<T^-1|H_x|T^0>, takes a minus sign and a transpose (because computed <T^0|H_x|T^-1>)
     747           8 :       CALL cp_fm_transpose(tmp_fm, work_fm)
     748           8 :       CALL cp_fm_scale(-1.0_dp, work_fm)
     749             :       CALL cp_fm_to_fm_submat(msource=work_fm, &
     750             :                               mtarget=img_fm, &
     751             :                               nrow=nex, &
     752             :                               ncol=nex, &
     753             :                               s_firstrow=1, &
     754             :                               s_firstcol=1, &
     755             :                               t_firstrow=1 + nsg + 1, &
     756           8 :                               t_firstcol=1 + nsg + ntp + 1)
     757             : 
     758             :       IF (debug_this_module) THEN
     759             :          WRITE (log_unit, "(/,A)") "<T^-1|H_x|T^0> in Hartree / cm^-1"
     760             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     761             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     762             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     763             :       END IF
     764             : 
     765             :       !Ms=0 to Ms=+-1 SOC, real part
     766             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     767             :                           orb_soc_y, dbcsr_tp, &
     768             :                           0.0_dp, dbcsr_work, &
     769           8 :                           filter_eps=eps_filter)
     770             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     771             :                           dbcsr_tp, dbcsr_work, &
     772             :                           0.0_dp, dbcsr_prod, &
     773           8 :                           filter_eps=eps_filter)
     774             : 
     775             :       CALL rcs_amew_soc_elements(dbcsr_tmp, &
     776             :                                  dbcsr_prod, &
     777             :                                  dbcsr_ovlp, &
     778             :                                  mo_soc_y, &
     779             :                                  pref_trace=1.0_dp, &
     780           8 :                                  pref_overall=0.5_dp*sqrt2)
     781             : 
     782             :       !<T^0|H_y|T^+1>
     783           8 :       CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
     784             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     785             :                               mtarget=real_fm, &
     786             :                               nrow=nex, &
     787             :                               ncol=nex, &
     788             :                               s_firstrow=1, &
     789             :                               s_firstcol=1, t_firstrow=1 + nsg + ntp + 1, &
     790           8 :                               t_firstcol=1 + nsg + 2*ntp + 1)
     791             : 
     792             :       IF (debug_this_module) THEN
     793             :          WRITE (log_unit, "(/,A)") "<T^0|H_y|T^+1> in Hartree / cm^-1"
     794             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     795             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     796             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     797             :       END IF
     798             : 
     799             :       !<T^-1|H_y|T^0>, takes a minus sign and a transpose
     800           8 :       CALL cp_fm_transpose(tmp_fm, work_fm)
     801           8 :       CALL cp_fm_scale(-1.0_dp, work_fm)
     802             :       CALL cp_fm_to_fm_submat(msource=work_fm, &
     803             :                               mtarget=real_fm, &
     804             :                               nrow=nex, &
     805             :                               ncol=nex, &
     806             :                               s_firstrow=1, &
     807             :                               s_firstcol=1, t_firstrow=1 + nsg + 1, &
     808           8 :                               t_firstcol=1 + nsg + ntp + 1)
     809             : 
     810             :       IF (debug_this_module) THEN
     811             :          WRITE (log_unit, "(/,A)") "<T^0|H_y|T^+1> in Hartree / cm^-1"
     812             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     813             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     814             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     815             :       END IF
     816             : 
     817             :       !Ms=1 to Ms=1 and Ms=-1 to Ms=-1 SOC, purely imaginary
     818             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, &
     819             :                           orb_soc_z, dbcsr_tp, &
     820             :                           0.0_dp, dbcsr_work, &
     821           8 :                           filter_eps=eps_filter)
     822             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, &
     823             :                           dbcsr_tp, dbcsr_work, &
     824             :                           0.0_dp, dbcsr_prod, &
     825           8 :                           filter_eps=eps_filter)
     826             : 
     827             :       CALL rcs_amew_soc_elements(dbcsr_tmp, &
     828             :                                  dbcsr_prod, &
     829             :                                  dbcsr_ovlp, &
     830             :                                  mo_soc_z, &
     831             :                                  pref_trace=1.0_dp, &
     832           8 :                                  pref_overall=1.0_dp)
     833             : 
     834             :       !<T^+1|H_z|T^+1>
     835           8 :       CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
     836             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     837             :                               mtarget=img_fm, &
     838             :                               nrow=nex, &
     839             :                               ncol=nex, &
     840             :                               s_firstrow=1, &
     841             :                               s_firstcol=1, &
     842             :                               t_firstrow=1 + nsg + 2*ntp + 1, &
     843           8 :                               t_firstcol=1 + nsg + 2*ntp + 1)
     844             : 
     845             :       IF (debug_this_module) THEN
     846             :          WRITE (log_unit, "(/,A)") "<T^+1|H_z|T^+1> in Hartree / cm^-1"
     847             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     848             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     849             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     850             :       END IF
     851             : 
     852             :       !<T^-1|H_z|T^-1>, takes a minus sign
     853           8 :       CALL cp_fm_scale(-1.0_dp, tmp_fm)
     854             :       CALL cp_fm_to_fm_submat(msource=tmp_fm, &
     855             :                               mtarget=img_fm, &
     856             :                               nrow=nex, &
     857             :                               ncol=nex, &
     858             :                               s_firstrow=1, &
     859             :                               s_firstcol=1, &
     860             :                               t_firstrow=1 + nsg + 1, &
     861           8 :                               t_firstcol=1 + nsg + 1)
     862             : 
     863             :       IF (debug_this_module) THEN
     864             :          WRITE (log_unit, "(/,A)") "<T^-1|H_z|T^-1> in Hartree / cm^-1"
     865             :          CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
     866             :          CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
     867             :          CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
     868             :       END IF
     869             : 
     870             :       !Intermediate clean-up
     871           8 :       CALL cp_fm_struct_release(work_struct)
     872           8 :       CALL cp_fm_release(work_fm)
     873           8 :       CALL cp_fm_release(tmp_fm)
     874           8 :       DEALLOCATE (diag, mo_soc_x, mo_soc_y, mo_soc_z)
     875             : 
     876           8 :       IF (log_unit > 0) THEN
     877           4 :          WRITE (log_unit, '(A)') "Diagonlzing SOC-Matrix"
     878             :       END IF
     879             : 
     880             :       !Set-up the complex hermitian perturbation matrix
     881           8 :       CALL cp_cfm_create(hami_cfm, full_struct)
     882           8 :       CALL cp_fm_to_cfm(real_fm, img_fm, hami_cfm)
     883             : 
     884             :       !!Optinal Output: SOC-Matrix
     885           8 :       IF (logger%para_env%is_source()) THEN
     886             :          log_unit = cp_print_key_unit_nr(logger, &
     887             :                                          tddfpt_print_section, &
     888             :                                          "SOC_PRINT", &
     889             :                                          extension=".socme", &
     890             :                                          file_form="FORMATTED", &
     891             :                                          file_action="WRITE", &
     892           4 :                                          file_status="UNKNOWN")
     893             :       ELSE
     894           4 :          log_unit = -1
     895             :       END IF
     896             : 
     897             :       !! Get the requested energy unit and optional outputs
     898           8 :       CALL section_vals_val_get(soc_print_section, "UNIT_eV", l_val=print_ev)
     899           8 :       CALL section_vals_val_get(soc_print_section, "UNIT_wn", l_val=print_wn)
     900           8 :       CALL section_vals_val_get(soc_print_section, "SPLITTING", l_val=print_splitting)
     901           8 :       CALL section_vals_val_get(soc_print_section, "SOME", l_val=print_some)
     902             : 
     903             :       !! save convertion unit into different varible for cleaner code
     904           8 :       IF (print_wn) THEN
     905           2 :          print_ev = .FALSE.
     906           2 :          unit_scale = wavenumbers
     907           6 :       ELSE IF (print_ev) THEN
     908             :          unit_scale = evolt
     909             :       ELSE
     910           0 :          CPWARN("No output unit was selected, printout will be in Hartree!")
     911           0 :          unit_scale = 1
     912             :       END IF
     913             : 
     914             :       !! This needs something to deal with a parallel run
     915             :       !! This output will be needed for NEWTONX for ISC!!
     916           8 :       IF (print_some) THEN
     917           2 :          IF (log_unit > 0) THEN
     918           1 :             WRITE (log_unit, '(A)') "____________________________________________________"
     919           1 :             WRITE (log_unit, '(A)') "                   FULL SOC-Matrix                  "
     920           1 :             IF (print_ev) THEN
     921           0 :                WRITE (log_unit, '(A)') "STATE   STATE   REAL-PART[eV]       IMG-PART[eV]    "
     922           1 :             ELSE IF (print_wn) THEN
     923           1 :                WRITE (log_unit, '(A)') "STATE   STATE   REAL-PART[cm-1]       IMG-PART[cm-1]"
     924             :             ELSE
     925           0 :                WRITE (log_unit, '(A)') "STATE   STATE   REAL-PART[Hartree] IMG-PART[Hartree]"
     926             :             END IF
     927           1 :             WRITE (log_unit, '(A)') "____________________________________________________"
     928             :          END IF
     929          12 :          ALLOCATE (real_tmp(ntot, ntot), img_tmp(ntot, ntot))
     930         182 :          real_tmp = 0_dp
     931         182 :          img_tmp = 0_dp
     932             : 
     933           2 :          CALL cp_fm_get_submatrix(real_fm, real_tmp, 1, 1, ntot, ntot)
     934           2 :          CALL cp_fm_get_submatrix(img_fm, img_tmp, 1, 1, ntot, ntot)
     935             : 
     936           2 :          IF (log_unit > 0) THEN
     937          10 :             DO jj = 1, ntot
     938          91 :                DO ii = 1, ntot
     939          81 :                   WRITE (unit=log_unit, fmt="(I3,5X,I3,5X,f16.8,5X,f16.8)") ii, jj, &
     940          81 :                      real_tmp(ii, jj)*unit_scale, &
     941         171 :                      img_tmp(ii, jj)*unit_scale
     942             :                END DO !! jj
     943             :             END DO        !! ii
     944           1 :             WRITE (log_unit, '(A)') "                   END SOC-MATRIX                   "
     945           1 :             WRITE (log_unit, '(A)') "____________________________________________________"
     946           1 :             WRITE (log_unit, '(A)') "____________________________________________________"
     947             :          END IF !(log_unit)
     948           2 :          DEALLOCATE (real_tmp, img_tmp)
     949             :       END IF !print_some
     950             : 
     951           8 :       CALL cp_fm_release(real_fm)
     952           8 :       CALL cp_fm_release(img_fm)
     953             : 
     954             :       ! Diagonalize the Hamiltonian
     955          24 :       ALLOCATE (tmp_evals(ntot))
     956           8 :       CALL cp_cfm_create(evecs_cfm, full_struct)
     957           8 :       CALL cp_cfm_heevd(hami_cfm, evecs_cfm, tmp_evals)
     958             : 
     959             :       !  Adjust the energies so the GS has zero, and store in the donor_state (without the GS)
     960          24 :       ALLOCATE (soc_env%soc_evals(ntot - 1))
     961          96 :       soc_env%soc_evals(:) = tmp_evals(2:ntot) - tmp_evals(1)
     962             : 
     963             :       !! We may be interested in the ground state stabilisation energy
     964           8 :       IF (logger%para_env%is_source()) THEN
     965           4 :          log_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     966             :       ELSE
     967           4 :          log_unit = -1
     968             :       END IF
     969             : 
     970           8 :       IF (log_unit > 0) THEN
     971           4 :          WRITE (log_unit, '(A27,6X,F18.10)') "Ground state stabilisation:", (soc_env%soc_evals(1)*unit_scale)
     972             :          IF (debug_this_module) THEN
     973             :             WRITE (log_unit, '(A)') "Calculation Dipole Moments"
     974             :          END IF
     975             :       END IF
     976             : 
     977             :       !Compute the dipole oscillator strengths
     978           8 :       soc_env%evals_a => evals_sing
     979           8 :       soc_env%evals_b => evals_trip
     980             :       CALL soc_dipol(qs_env, dbcsr_soc_package, soc_env, &
     981             :                      evecs_cfm, eps_filter, dipole_form, gs_mos, &
     982           8 :                      evects_sing, evects_trip)
     983             : 
     984             :       !Create final output
     985             :       !! Output unit is choosen by the keyword "UNIT_eV" and "UNIT_wn" in "SOC_PRINT" section
     986           8 :       IF (logger%para_env%is_source()) THEN
     987           4 :          log_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     988             :       ELSE
     989           4 :          log_unit = -1
     990             :       END IF
     991             : 
     992           8 :       IF (log_unit > 0) THEN
     993           4 :          WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
     994           4 :          WRITE (log_unit, '(A)') "-                         SOC CORRECTED SPECTRUM                             -"
     995           4 :          WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
     996           4 :          IF (print_ev) THEN
     997           3 :             WRITE (log_unit, '(A)') "-     STATE   SOC-corrected exc. energies [eV]  Oscillator strengths [a.u.]  -"
     998           1 :          ELSE IF (print_wn) THEN
     999           1 :             WRITE (log_unit, '(A)') "-     STATE   SOC-corrected exc. energies [cm-1] Oscillator strengths [a.u.]-"
    1000             :          ELSE
    1001           0 :             WRITE (log_unit, '(A)') "-     STATE SOC-corrected exc.energies [Hartree] Oscillator strengths [a.u.]-"
    1002             :          END IF
    1003           4 :          WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
    1004          48 :          DO istate = 1, SIZE(soc_env%soc_evals)
    1005          44 :             WRITE (log_unit, '(6X,I5,11X,2F16.5)') istate, soc_env%soc_evals(istate)*unit_scale, &
    1006          92 :                soc_env%soc_osc(istate)
    1007             :          END DO
    1008             :          !! This is used for regtesting
    1009          48 :          WRITE (log_unit, '(A,F18.8)') "TDDFT+SOC : CheckSum (eV) ", SUM(soc_env%soc_evals)*evolt
    1010             :       END IF
    1011             : 
    1012             :       !! We may be interested in the differenz between the SOC and
    1013             :       !! TDDFPT Eigenvalues
    1014             :       !! Activated with the "SPLITTING" keyword in "SOC_PRINT" section
    1015          16 :       IF (print_splitting .OR. debug_this_module) THEN
    1016           8 :          CALL resort_evects(evecs_cfm, evects_sort)
    1017           8 :          IF (log_unit > 0) THEN
    1018           4 :             WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
    1019           4 :             WRITE (log_unit, '(A)') "-                              Splittings                                   -"
    1020           4 :             WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
    1021           4 :             IF (print_ev) THEN
    1022           3 :                WRITE (log_unit, '(A)') "-     STATE          exc.energies[eV]          splittings[eV]               -"
    1023           1 :             ELSE IF (print_wn) THEN
    1024           1 :                WRITE (log_unit, '(A)') "-     STATE          exc.energies[cm-1]        splittings[cm-1]             -"
    1025             :             ELSE
    1026           0 :                WRITE (log_unit, '(A)') "-     STATE          exc.energies[Hartree]     splittings[Hartree]          -"
    1027             :             END IF
    1028           4 :             WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
    1029           4 :             mult = "   "
    1030          48 :             DO iex = 1, SIZE(soc_env%soc_evals)
    1031          44 :                IF (evects_sort(iex + 1) .GT. nsg + ntp*2 + 1) THEN
    1032          11 :                   mult = "T+1"
    1033          33 :                ELSE IF (evects_sort(iex + 1) .GT. nsg + ntp + 1) THEN
    1034          11 :                   mult = "T0"
    1035          22 :                ELSE IF (evects_sort(iex + 1) .GT. nsg + 1) THEN
    1036          11 :                   mult = "T-1"
    1037             :                ELSE
    1038          11 :                   mult = "S  "
    1039             :                END IF
    1040          48 :                IF (mult == "T-1") THEN
    1041          11 :                   WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
    1042          11 :                      evects_sort(iex + 1) - (1 + nsg), soc_env%soc_evals(iex)*unit_scale, &
    1043          22 :                      (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg)))*unit_scale
    1044          33 :                ELSE IF (mult == "T0") THEN
    1045          11 :                   WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
    1046          11 :                      evects_sort(iex + 1) - (1 + nsg + ntp), soc_env%soc_evals(iex)*unit_scale, &
    1047          22 :                      (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg + ntp)))*unit_scale
    1048          22 :                ELSE IF (mult == "T+1") THEN
    1049          11 :                   WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
    1050          11 :                      evects_sort(iex + 1) - (1 + nsg + ntp*2), soc_env%soc_evals(iex)*unit_scale, &
    1051          22 :                      (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg + ntp*2)))*unit_scale
    1052          11 :                ELSE IF (mult == "S  ") THEN
    1053          11 :                   WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
    1054          11 :                      evects_sort(iex + 1), soc_env%soc_evals(iex)*unit_scale, &
    1055          22 :                      (soc_env%soc_evals(iex) - evals_sing(evects_sort(iex + 1) - 1))*unit_scale
    1056             :                END IF
    1057             :             END DO
    1058           4 :             WRITE (log_unit, '(A)') "----------------------------------------------------------------------------"
    1059           4 :             DEALLOCATE (evects_sort)
    1060             :          END IF
    1061             :       END IF
    1062             : 
    1063             :       !! clean up
    1064           8 :       CALL soc_env_release(soc_env)
    1065           8 :       CALL cp_fm_struct_release(full_struct)
    1066           8 :       CALL cp_cfm_release(hami_cfm)
    1067           8 :       CALL cp_cfm_release(evecs_cfm)
    1068           8 :       CALL dbcsr_distribution_release(coeffs_dist)
    1069           8 :       CALL dbcsr_distribution_release(prod_dist)
    1070           8 :       CALL dbcsr_release(dbcsr_sg)
    1071           8 :       CALL dbcsr_release(dbcsr_tp)
    1072           8 :       CALL dbcsr_release(dbcsr_prod)
    1073           8 :       CALL dbcsr_release(dbcsr_ovlp)
    1074           8 :       CALL dbcsr_release(dbcsr_tmp)
    1075           8 :       CALL dbcsr_release(dbcsr_work)
    1076             :       IF (debug_this_module) THEN
    1077             :          CALL dbcsr_release(dbcsr_dummy)
    1078             :          DEALLOCATE (dbcsr_dummy)
    1079             :       END IF
    1080             : 
    1081           8 :       DEALLOCATE (dbcsr_work, dbcsr_prod, dbcsr_ovlp, dbcsr_tmp)
    1082           8 :       DEALLOCATE (coeffs_dist, prod_dist, col_dist, col_blk_size, row_dist_new)
    1083           8 :       DEALLOCATE (dbcsr_sg, dbcsr_tp, tmp_evals)
    1084             : 
    1085           8 :       CALL timestop(handle)
    1086             : 
    1087         112 :    END SUBROUTINE tddfpt_soc_rcs
    1088             : 
    1089             : ! **************************************************************************************************
    1090             : !> \brief Some additional initalizations
    1091             : !> \param qs_env Quickstep environment
    1092             : !> \param soc_atom_env contains information to build the ZORA-Hamiltionian
    1093             : !> \param soc_env contails details and outputs for SOC Correction
    1094             : !> \param evects_a singlet Eigenvectors
    1095             : !> \param evects_b triplet eigenvectors
    1096             : !> \param dipole_form choosen dipole-form from TDDFPT-Section
    1097             : !> \param gs_mos ...
    1098             : ! **************************************************************************************************
    1099           8 :    SUBROUTINE inititialize_soc(qs_env, soc_atom_env, soc_env, &
    1100           8 :                                evects_a, evects_b, &
    1101           8 :                                dipole_form, gs_mos)
    1102             :       !! Different Types
    1103             :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
    1104             :       TYPE(soc_atom_env_type), INTENT(OUT), POINTER      :: soc_atom_env
    1105             :       TYPE(soc_env_type), INTENT(OUT), TARGET            :: soc_env
    1106             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects_a, evects_b
    1107             :       INTEGER                                            :: dipole_form
    1108             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1109             :          INTENT(in)                                      :: gs_mos
    1110             : 
    1111             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'inititialize_soc'
    1112             : 
    1113             :       CHARACTER(len=default_string_length), &
    1114             :          ALLOCATABLE, DIMENSION(:, :)                    :: grid_info
    1115             :       CHARACTER(len=default_string_length), &
    1116           8 :          DIMENSION(:), POINTER                           :: k_list
    1117             :       INTEGER                                            :: handle, i, ikind, irep, ispin, nactive, &
    1118             :                                                             nao, natom, nkind, nrep, nspins, &
    1119             :                                                             nstates
    1120             :       LOGICAL                                            :: all_potential_present, do_os
    1121             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1122             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
    1123             :       TYPE(cp_fm_type), POINTER                          :: a_coeff, b_coeff
    1124           8 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
    1125             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1126           8 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1127             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1128           8 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1129           8 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1130             :       TYPE(section_vals_type), POINTER                   :: soc_section
    1131             :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
    1132             : 
    1133           8 :       CALL timeset(routineN, handle)
    1134             : 
    1135           8 :       NULLIFY (qs_kind_set, particle_set, blacs_env, para_env, &
    1136           8 :                a_coeff, b_coeff, tddfpt_control, soc_section)
    1137             : 
    1138             :       !! Open_shell is not implemented yet
    1139           8 :       do_os = .FALSE.
    1140             : 
    1141             :       !! soc_env will pass most of the larger matrices to the different modules
    1142           8 :       CALL soc_env_create(soc_env)
    1143             :       CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, &
    1144             :                       mos=mos, natom=natom, &
    1145             :                       particle_set=particle_set, blacs_env=blacs_env, &
    1146             :                       para_env=para_env, matrix_s=matrix_s, &
    1147           8 :                       dft_control=dft_control)
    1148             : 
    1149             :       !! The Dipole form shall be consistent with the tddfpt-module!
    1150           8 :       tddfpt_control => dft_control%tddfpt2_control
    1151           8 :       dipole_form = tddfpt_control%dipole_form
    1152             : 
    1153           8 :       soc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%SOC")
    1154             : 
    1155             :       !! We may want to use a bigger grid then default
    1156           8 :       CALL section_vals_val_get(soc_section, "GRID", n_rep_val=nrep)
    1157          16 :       ALLOCATE (grid_info(nrep, 3))
    1158           8 :       DO irep = 1, nrep
    1159             :          CALL section_vals_val_get(soc_section, &
    1160             :                                    "GRID", &
    1161             :                                    i_rep_val=irep, &
    1162           0 :                                    c_vals=k_list)
    1163           8 :          grid_info(irep, :) = k_list
    1164             :       END DO
    1165             : 
    1166           8 :       CALL soc_dipole_operator(soc_env, tddfpt_control, qs_env, gs_mos)
    1167             : 
    1168           8 :       nspins = SIZE(evects_a, 1)
    1169          16 :       DO ispin = 1, nspins
    1170          16 :          CALL cp_fm_get_info(evects_a(ispin, 1), nrow_global=nao, ncol_global=nactive)
    1171             :       END DO
    1172             : 
    1173             :       !! We need to change and create structures for the exitation vectors.
    1174             :       !! All excited states shall be written in one matrix, oneafter the other
    1175           8 :       nstates = SIZE(evects_a, 2)
    1176             : 
    1177           8 :       a_coeff => soc_env%a_coeff
    1178           8 :       b_coeff => soc_env%b_coeff
    1179             : 
    1180           8 :       NULLIFY (fm_struct)
    1181             :       CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
    1182           8 :                                nrow_global=nao, ncol_global=nspins*nactive*nstates)
    1183           8 :       CALL cp_fm_create(a_coeff, fm_struct)
    1184           8 :       CALL cp_fm_create(b_coeff, fm_struct)
    1185           8 :       CALL cp_fm_struct_release(fm_struct)
    1186             : 
    1187           8 :       CALL soc_contract_evect(evects_a, a_coeff)
    1188           8 :       CALL soc_contract_evect(evects_b, b_coeff)
    1189             : 
    1190          32 :       ALLOCATE (soc_env%orb_soc(3))
    1191          32 :       DO i = 1, 3
    1192          32 :          ALLOCATE (soc_env%orb_soc(i)%matrix)
    1193             :       END DO
    1194             : 
    1195             :       ! Make soc_atom_env for H^SOC calculations
    1196             :       ! Compute the contraction coefficients for spherical orbitals
    1197           8 :       CALL soc_atom_create(soc_atom_env)
    1198             :       IF (do_os) THEN
    1199             :          soc_atom_env%nspins = 2
    1200             :       ELSE
    1201           8 :          soc_atom_env%nspins = 1
    1202             :       END IF
    1203             : 
    1204           8 :       nkind = SIZE(qs_kind_set)
    1205          32 :       ALLOCATE (soc_atom_env%grid_atom_set(nkind))
    1206          32 :       ALLOCATE (soc_atom_env%harmonics_atom_set(nkind))
    1207          32 :       ALLOCATE (soc_atom_env%orb_sphi_so(nkind))
    1208             : 
    1209           8 :       CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
    1210           8 :       IF (.NOT. all_potential_present) THEN
    1211           2 :          CALL V_SOC_xyz_from_pseudopotential(qs_env, soc_atom_env%soc_pp)
    1212             :       END IF
    1213             : 
    1214             :       !! Prepare the atomic grid we need to calculate the SOC Operator in an Atomic basis
    1215             :       !! copied from init_xas_atom_grid_harmo for convenience
    1216           8 :       CALL init_atom_grid(soc_atom_env, grid_info, qs_env)
    1217             : 
    1218          16 :       DO ikind = 1, nkind
    1219           8 :          NULLIFY (soc_atom_env%orb_sphi_so(ikind)%array)
    1220          16 :          CALL compute_sphi_so(ikind, "ORB", soc_atom_env%orb_sphi_so(ikind)%array, qs_env)
    1221             :       END DO !ikind
    1222             : 
    1223           8 :       DEALLOCATE (grid_info)
    1224             : 
    1225           8 :       CALL timestop(handle)
    1226             : 
    1227          40 :    END SUBROUTINE inititialize_soc
    1228             : 
    1229             : ! **************************************************************************************************
    1230             : !> \brief Initializes the atomic grids and harmonics for the atomic calculations
    1231             : !> \param soc_atom_env ...
    1232             : !> \param grid_info ...
    1233             : !> \param qs_env ...
    1234             : !> \note Copied and modified from init_xas_atom_grid_harmo
    1235             : ! **************************************************************************************************
    1236           8 :    SUBROUTINE init_atom_grid(soc_atom_env, grid_info, qs_env)
    1237             : 
    1238             :       TYPE(soc_atom_env_type)                            :: soc_atom_env
    1239             :       CHARACTER(len=default_string_length), &
    1240             :          ALLOCATABLE, DIMENSION(:, :)                    :: grid_info
    1241             :       TYPE(qs_environment_type)                          :: qs_env
    1242             : 
    1243             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'init_atom_grid'
    1244             : 
    1245             :       CHARACTER(LEN=default_string_length)               :: kind_name
    1246             :       INTEGER :: handle, igrid, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, &
    1247             :          llmax, lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nr, &
    1248             :          quadrature, stat
    1249             :       REAL(dp)                                           :: kind_radius
    1250             :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: rga
    1251           8 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: my_CG
    1252             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1253             :       TYPE(grid_atom_type), POINTER                      :: grid_atom
    1254             :       TYPE(gto_basis_set_type), POINTER                  :: tmp_basis
    1255             :       TYPE(harmonics_atom_type), POINTER                 :: harmonics
    1256             :       TYPE(qs_control_type), POINTER                     :: qs_control
    1257           8 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1258             : 
    1259           8 :       CALL timeset(routineN, handle)
    1260             : 
    1261           8 :       NULLIFY (my_CG, qs_kind_set, dft_control, qs_control)
    1262           8 :       NULLIFY (grid_atom, harmonics, tmp_basis)
    1263             : 
    1264             :       !!  Initialization of some integer for the CG coeff generation
    1265             :       CALL get_qs_env(qs_env, &
    1266             :                       qs_kind_set=qs_kind_set, &
    1267           8 :                       dft_control=dft_control)
    1268             : 
    1269           8 :       CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="ORB")
    1270           8 :       qs_control => dft_control%qs_control
    1271             : 
    1272             :       !!To be sure:: Is the basis grid present?
    1273           8 :       IF (maxlgto < 0) THEN
    1274           0 :          CPABORT("tmp_basis is not Present!")
    1275             :       END IF
    1276             : 
    1277             :       !maximum expansion
    1278           8 :       llmax = 2*maxlgto
    1279           8 :       max_s_harm = nsoset(llmax)
    1280           8 :       max_s_set = nsoset(maxlgto)
    1281           8 :       lcleb = llmax
    1282             : 
    1283             :       !  Allocate and compute the CG coeffs (copied from init_rho_atom)
    1284           8 :       CALL clebsch_gordon_init(lcleb)
    1285           8 :       CALL reallocate(my_CG, 1, max_s_set, 1, max_s_set, 1, max_s_harm)
    1286             : 
    1287          24 :       ALLOCATE (rga(lcleb, 2))
    1288          32 :       DO lc1 = 0, maxlgto
    1289         104 :          DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
    1290          72 :             l1 = indso(1, iso1)
    1291          72 :             m1 = indso(2, iso1)
    1292         312 :             DO lc2 = 0, maxlgto
    1293         936 :                DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
    1294         648 :                   l2 = indso(1, iso2)
    1295         648 :                   m2 = indso(2, iso2)
    1296         648 :                   CALL clebsch_gordon(l1, m1, l2, m2, rga)
    1297         648 :                   IF (l1 + l2 > llmax) THEN
    1298             :                      l1l2 = llmax
    1299             :                   ELSE
    1300             :                      l1l2 = l1 + l2
    1301             :                   END IF
    1302         648 :                   mp = m1 + m2
    1303         648 :                   mm = m1 - m2
    1304         648 :                   IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
    1305         288 :                      mp = -ABS(mp)
    1306         288 :                      mm = -ABS(mm)
    1307             :                   ELSE
    1308         360 :                      mp = ABS(mp)
    1309         360 :                      mm = ABS(mm)
    1310             :                   END IF
    1311        2304 :                   DO lp = MOD(l1 + l2, 2), l1l2, 2
    1312        1440 :                      il = lp/2 + 1
    1313        1440 :                      IF (ABS(mp) <= lp) THEN
    1314        1024 :                      IF (mp >= 0) THEN
    1315         688 :                         iso = nsoset(lp - 1) + lp + 1 + mp
    1316             :                      ELSE
    1317         336 :                         iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
    1318             :                      END IF
    1319        1024 :                      my_CG(iso1, iso2, iso) = rga(il, 1)
    1320             :                      END IF
    1321        2088 :                      IF (mp /= mm .AND. ABS(mm) <= lp) THEN
    1322         480 :                      IF (mm >= 0) THEN
    1323         320 :                         iso = nsoset(lp - 1) + lp + 1 + mm
    1324             :                      ELSE
    1325         160 :                         iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
    1326             :                      END IF
    1327         480 :                      my_CG(iso1, iso2, iso) = rga(il, 2)
    1328             :                      END IF
    1329             :                   END DO
    1330             :                END DO ! iso2
    1331             :             END DO ! lc2
    1332             :          END DO ! iso1
    1333             :       END DO ! lc1
    1334           8 :       DEALLOCATE (rga)
    1335           8 :       CALL clebsch_gordon_deallocate()
    1336             : 
    1337             :       !  Create the Lebedev grids and compute the spherical harmonics
    1338           8 :       CALL init_lebedev_grids()
    1339           8 :       quadrature = qs_control%gapw_control%quadrature
    1340             : 
    1341          16 :       DO ikind = 1, SIZE(soc_atom_env%grid_atom_set)
    1342             : 
    1343             :          ! Allocate the grid and the harmonics for this kind
    1344           8 :          NULLIFY (soc_atom_env%grid_atom_set(ikind)%grid_atom)
    1345           8 :          NULLIFY (soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
    1346           8 :          CALL allocate_grid_atom(soc_atom_env%grid_atom_set(ikind)%grid_atom)
    1347             :          CALL allocate_harmonics_atom( &
    1348           8 :             soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
    1349             : 
    1350             :          NULLIFY (grid_atom, harmonics)
    1351           8 :          grid_atom => soc_atom_env%grid_atom_set(ikind)%grid_atom
    1352           8 :          harmonics => soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom
    1353             : 
    1354             :          ! Initialize some integers
    1355             :          CALL get_qs_kind(qs_kind_set(ikind), &
    1356             :                           ngrid_rad=nr, &
    1357             :                           ngrid_ang=na, &
    1358           8 :                           name=kind_name)
    1359             : 
    1360             :          ! take the grid dimension given as input, if none, take the GAPW ones above
    1361           8 :          DO igrid = 1, SIZE(grid_info, 1)
    1362           8 :             IF (grid_info(igrid, 1) == kind_name) THEN
    1363           0 :                READ (grid_info(igrid, 2), *, iostat=stat) na
    1364           0 :                IF (stat .NE. 0) CPABORT("The 'na' value for the GRID keyword must be an integer")
    1365           0 :                READ (grid_info(igrid, 3), *, iostat=stat) nr
    1366           0 :                IF (stat .NE. 0) CPABORT("The 'nr' value for the GRID keyword must be an integer")
    1367             :                EXIT
    1368             :             END IF
    1369             :          END DO !igrid
    1370             : 
    1371           8 :          ll = get_number_of_lebedev_grid(n=na)
    1372           8 :          na = lebedev_grid(ll)%n
    1373           8 :          la = lebedev_grid(ll)%l
    1374           8 :          grid_atom%ng_sphere = na
    1375           8 :          grid_atom%nr = nr
    1376             : 
    1377             :          !Create the harmonics with the ORB-Basis
    1378             :          CALL get_qs_kind(qs_kind_set(ikind), &
    1379             :                           basis_set=tmp_basis, &
    1380           8 :                           basis_type="ORB")
    1381             :          CALL get_gto_basis_set(gto_basis_set=tmp_basis, &
    1382             :                                 maxl=maxl, &
    1383           8 :                                 kind_radius=kind_radius)
    1384             : 
    1385           8 :          CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
    1386           8 :          CALL truncate_radial_grid(grid_atom, kind_radius)
    1387             : 
    1388           8 :          maxs = nsoset(maxl)
    1389             :          CALL create_harmonics_atom(harmonics, &
    1390             :                                     my_CG, na, &
    1391             :                                     llmax, maxs, &
    1392             :                                     max_s_harm, ll, &
    1393             :                                     grid_atom%wa, &
    1394             :                                     grid_atom%azi, &
    1395           8 :                                     grid_atom%pol)
    1396          32 :          CALL get_maxl_CG(harmonics, tmp_basis, llmax, max_s_harm)
    1397             :       END DO !ikind
    1398             : 
    1399           8 :       CALL deallocate_lebedev_grids()
    1400           8 :       DEALLOCATE (my_CG)
    1401             : 
    1402           8 :       CALL timestop(handle)
    1403             : 
    1404          16 :    END SUBROUTINE init_atom_grid
    1405             : 
    1406             : ! **************************************************************************************************
    1407             : !> \brief ...
    1408             : !> \param qs_env ...
    1409             : !> \param dbcsr_soc_package ...
    1410             : !> \param soc_env ...
    1411             : !> \param soc_evecs_cfm ...
    1412             : !> \param eps_filter ...
    1413             : !> \param dipole_form ...
    1414             : !> \param gs_mos ...
    1415             : !> \param evects_sing ...
    1416             : !> \param evects_trip ...
    1417             : ! **************************************************************************************************
    1418           8 :    SUBROUTINE soc_dipol(qs_env, dbcsr_soc_package, soc_env, &
    1419             :                         soc_evecs_cfm, eps_filter, dipole_form, gs_mos, &
    1420           8 :                         evects_sing, evects_trip)
    1421             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1422             :       TYPE(dbcsr_soc_package_type)                       :: dbcsr_soc_package
    1423             :       TYPE(soc_env_type), TARGET                         :: soc_env
    1424             :       TYPE(cp_cfm_type), INTENT(in)                      :: soc_evecs_cfm
    1425             :       REAL(dp), INTENT(IN)                               :: eps_filter
    1426             :       INTEGER                                            :: dipole_form
    1427             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1428             :          POINTER                                         :: gs_mos
    1429             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects_sing, evects_trip
    1430             : 
    1431             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'soc_dipol'
    1432             : 
    1433             :       COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :)          :: transdip
    1434             :       INTEGER                                            :: handle, i, nosc, ntot
    1435           8 :       REAL(dp), DIMENSION(:), POINTER                    :: osc_str, soc_evals
    1436             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1437             :       TYPE(cp_cfm_type)                                  :: dip_cfm, work1_cfm, work2_cfm
    1438             :       TYPE(cp_fm_struct_type), POINTER                   :: dip_struct, full_struct
    1439           8 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: amew_dip
    1440             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1441             : 
    1442           8 :       CALL timeset(routineN, handle)
    1443             : 
    1444           8 :       NULLIFY (para_env, blacs_env, dip_struct, full_struct, osc_str)
    1445           8 :       NULLIFY (soc_evals)
    1446             : 
    1447           8 :       CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
    1448             : 
    1449           8 :       soc_evals => soc_env%soc_evals
    1450           8 :       nosc = SIZE(soc_evals)
    1451           8 :       ntot = nosc + 1
    1452          24 :       ALLOCATE (soc_env%soc_osc(nosc))
    1453           8 :       osc_str => soc_env%soc_osc
    1454          96 :       osc_str(:) = 0.0_dp
    1455             : 
    1456             :       !get some work arrays/matrix
    1457             :       CALL cp_fm_struct_create(dip_struct, &
    1458             :                                context=blacs_env, para_env=para_env, &
    1459           8 :                                nrow_global=ntot, ncol_global=1)
    1460             : 
    1461           8 :       CALL cp_cfm_get_info(soc_evecs_cfm, matrix_struct=full_struct)
    1462           8 :       CALL cp_cfm_create(dip_cfm, dip_struct)
    1463           8 :       CALL cp_cfm_create(work1_cfm, full_struct)
    1464           8 :       CALL cp_cfm_create(work2_cfm, full_struct)
    1465             : 
    1466          24 :       ALLOCATE (transdip(ntot, 1))
    1467             : 
    1468             :       CALL get_rcs_amew_op(amew_dip, soc_env, dbcsr_soc_package, &
    1469             :                            eps_filter, qs_env, gs_mos, &
    1470           8 :                            evects_sing, evects_trip)
    1471             : 
    1472          32 :       DO i = 1, 3 !cartesian coord x, y, z
    1473             : 
    1474             :          !Convert the real dipole into the cfm format for calculations
    1475          24 :          CALL cp_fm_to_cfm(msourcer=amew_dip(i), mtarget=work1_cfm)
    1476             : 
    1477             :          !compute amew_coeffs^dagger * amew_dip * amew_gs to get the transition moments
    1478             :          CALL parallel_gemm('C', 'N', ntot, &
    1479             :                             ntot, ntot, &
    1480             :                             (1.0_dp, 0.0_dp), &
    1481             :                             soc_evecs_cfm, &
    1482             :                             work1_cfm, &
    1483             :                             (0.0_dp, 0.0_dp), &
    1484          24 :                             work2_cfm)
    1485             :          CALL parallel_gemm('N', 'N', ntot, &
    1486             :                             1, ntot, &
    1487             :                             (1.0_dp, 0.0_dp), &
    1488             :                             work2_cfm, &
    1489             :                             soc_evecs_cfm, &
    1490             :                             (0.0_dp, 0.0_dp), &
    1491          24 :                             dip_cfm)
    1492             : 
    1493          24 :          CALL cp_cfm_get_submatrix(dip_cfm, transdip)
    1494             : 
    1495             :          !transition dipoles are real numbers
    1496         296 :          osc_str(:) = osc_str(:) + REAL(transdip(2:ntot, 1))**2 + AIMAG(transdip(2:ntot, 1))**2
    1497             : 
    1498             :       END DO !i
    1499             : 
    1500             :       !multiply with appropriate prefac depending in the rep
    1501           8 :       IF (dipole_form == tddfpt_dipole_length) THEN
    1502         156 :          osc_str(:) = 2.0_dp/3.0_dp*soc_evals(:)*osc_str(:)
    1503             :       ELSE
    1504          36 :          osc_str(:) = 2.0_dp/3.0_dp*soc_evals(:)*osc_str(:)
    1505             :       END IF
    1506             : 
    1507             :       !clean-up
    1508           8 :       CALL cp_fm_struct_release(dip_struct)
    1509           8 :       CALL cp_cfm_release(work1_cfm)
    1510           8 :       CALL cp_cfm_release(work2_cfm)
    1511           8 :       CALL cp_cfm_release(dip_cfm)
    1512          32 :       DO i = 1, 3
    1513          32 :          CALL cp_fm_release(amew_dip(i))
    1514             :       END DO
    1515           8 :       DEALLOCATE (amew_dip, transdip)
    1516             : 
    1517           8 :       CALL timestop(handle)
    1518             : 
    1519          16 :    END SUBROUTINE soc_dipol
    1520             : 
    1521             : !**********************************************************************************
    1522             : !> \brief: This will create the Dipole operator within the amew basis
    1523             : !> \param amew_op Output dipole operator
    1524             : !> \param soc_env ...
    1525             : !> \param dbcsr_soc_package Some work matrices
    1526             : !> \param eps_filter ...
    1527             : !> \param qs_env ...
    1528             : !> \param gs_mos ...
    1529             : !> \param evects_sing ...
    1530             : !> \param evects_trip ...
    1531             : ! **************************************************************************************************
    1532           8 :    SUBROUTINE get_rcs_amew_op(amew_op, soc_env, dbcsr_soc_package, eps_filter, qs_env, gs_mos, &
    1533           8 :                               evects_sing, evects_trip)
    1534             : 
    1535             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
    1536             :          INTENT(OUT)                                     :: amew_op
    1537             :       TYPE(soc_env_type), TARGET                         :: soc_env
    1538             :       TYPE(dbcsr_soc_package_type)                       :: dbcsr_soc_package
    1539             :       REAL(dp), INTENT(IN)                               :: eps_filter
    1540             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1541             :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1542             :          POINTER                                         :: gs_mos
    1543             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evects_sing, evects_trip
    1544             : 
    1545             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'get_rcs_amew_op'
    1546             : 
    1547             :       INTEGER                                            :: dim_op, handle, i, isg, nactive, nao, &
    1548             :                                                             nex, nsg, ntot, ntp
    1549             :       LOGICAL                                            :: vel_rep
    1550             :       REAL(dp)                                           :: op, sqrt2
    1551           8 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: diag, gs_diag, gsgs_op
    1552           8 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: mo_op, sggs_block
    1553             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1554             :       TYPE(cp_fm_struct_type), POINTER                   :: full_struct, gsgs_struct, sggs_struct, &
    1555             :                                                             std_struct, tmp_struct
    1556             :       TYPE(cp_fm_type)                                   :: gs_fm, sggs_fm, tmp_fm, vec_op, work_fm
    1557             :       TYPE(cp_fm_type), POINTER                          :: gs_coeffs
    1558           8 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ao_op, matrix_s
    1559             :       TYPE(dbcsr_type), POINTER                          :: ao_op_i, dbcsr_ovlp, dbcsr_prod, &
    1560             :                                                             dbcsr_sg, dbcsr_tmp, dbcsr_tp, &
    1561             :                                                             dbcsr_work
    1562           8 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1563             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1564             : 
    1565           8 :       CALL timeset(routineN, handle)
    1566             : 
    1567           8 :       NULLIFY (mos, gs_coeffs)
    1568           8 :       NULLIFY (matrix_s)
    1569           8 :       NULLIFY (para_env, blacs_env)
    1570           8 :       NULLIFY (full_struct, gsgs_struct, std_struct, tmp_struct, sggs_struct)
    1571           8 :       NULLIFY (ao_op_i, ao_op)
    1572           8 :       NULLIFY (dbcsr_tp, dbcsr_sg, dbcsr_ovlp, dbcsr_work, dbcsr_tmp, dbcsr_prod)
    1573             : 
    1574             :       ! Initialization
    1575             :       !! Fist case for length, secound for velocity-representation
    1576           8 :       IF (ASSOCIATED(soc_env%dipmat)) THEN
    1577           6 :          ao_op => soc_env%dipmat
    1578           6 :          vel_rep = .FALSE.
    1579             :       ELSE
    1580           2 :          ao_op => soc_env%dipmat_ao
    1581           2 :          vel_rep = .TRUE.
    1582             :       END IF
    1583           8 :       nsg = SIZE(soc_env%evals_a)
    1584           8 :       ntp = nsg; nex = nsg !all the same by construction, keep them separate for clarity
    1585           8 :       ntot = 1 + nsg + 3*ntp
    1586             : 
    1587           8 :       CALL get_qs_env(qs_env, matrix_s=matrix_s, mos=mos, para_env=para_env, blacs_env=blacs_env)
    1588           8 :       sqrt2 = SQRT(2.0_dp)
    1589           8 :       dim_op = SIZE(ao_op)
    1590             : 
    1591           8 :       dbcsr_sg => dbcsr_soc_package%dbcsr_sg
    1592           8 :       dbcsr_tp => dbcsr_soc_package%dbcsr_tp
    1593           8 :       dbcsr_work => dbcsr_soc_package%dbcsr_work
    1594           8 :       dbcsr_prod => dbcsr_soc_package%dbcsr_prod
    1595           8 :       dbcsr_ovlp => dbcsr_soc_package%dbcsr_ovlp
    1596           8 :       dbcsr_tmp => dbcsr_soc_package%dbcsr_tmp
    1597             : 
    1598             :       !  Create the amew_op matrix
    1599             :       CALL cp_fm_struct_create(full_struct, &
    1600             :                                context=blacs_env, para_env=para_env, &
    1601           8 :                                nrow_global=ntot, ncol_global=ntot)
    1602          48 :       ALLOCATE (amew_op(dim_op))
    1603          32 :       DO i = 1, dim_op
    1604          32 :          CALL cp_fm_create(amew_op(i), full_struct)
    1605             :       END DO !i
    1606             : 
    1607             :       !  Deal with the GS-GS contribution <0|0> = 2*sum_j <phi_j|op|phi_j>
    1608           8 :       CALL get_mo_set(mos(1), nao=nao, homo=nactive)
    1609           8 :       gs_coeffs => gs_mos(1)%mos_occ
    1610           8 :       CALL cp_fm_get_info(gs_coeffs, matrix_struct=std_struct)
    1611             : 
    1612             :       CALL cp_fm_struct_create(gsgs_struct, &
    1613             :                                context=blacs_env, para_env=para_env, &
    1614           8 :                                nrow_global=nactive, ncol_global=nactive)
    1615             : 
    1616           8 :       CALL cp_fm_create(gs_fm, gsgs_struct) !nactive nactive
    1617           8 :       CALL cp_fm_create(work_fm, std_struct) !nao nactive
    1618             : 
    1619          24 :       ALLOCATE (gsgs_op(dim_op))
    1620          24 :       ALLOCATE (gs_diag(nactive))
    1621             : 
    1622          32 :       DO i = 1, dim_op
    1623             : 
    1624          24 :          ao_op_i => ao_op(i)%matrix
    1625          24 :          CALL cp_dbcsr_sm_fm_multiply(ao_op_i, gs_coeffs, work_fm, ncol=nactive)
    1626             :          CALL parallel_gemm('T', 'N', nactive, nactive, nao, &
    1627          24 :                             1.0_dp, gs_coeffs, work_fm, 0.0_dp, gs_fm)
    1628          24 :          CALL cp_fm_get_diag(gs_fm, gs_diag)
    1629         194 :          gsgs_op(i) = 2.0_dp*SUM(gs_diag)
    1630             : 
    1631             :       END DO !i
    1632           8 :       DEALLOCATE (gs_diag)
    1633             : 
    1634           8 :       CALL cp_fm_release(work_fm)
    1635             : 
    1636             :       !  Create the work and helper fms
    1637           8 :       CALL cp_fm_create(vec_op, std_struct) !nao nactive
    1638           8 :       CALL cp_fm_struct_release(gsgs_struct)
    1639             : 
    1640             :       CALL cp_fm_struct_create(tmp_struct, &
    1641             :                                context=blacs_env, para_env=para_env, &
    1642           8 :                                nrow_global=nex, ncol_global=nex)
    1643             :       CALL cp_fm_struct_create(sggs_struct, &
    1644             :                                context=blacs_env, para_env=para_env, &
    1645           8 :                                nrow_global=nactive*nsg, ncol_global=nactive)
    1646             : 
    1647           8 :       CALL cp_fm_create(tmp_fm, tmp_struct)
    1648           8 :       CALL cp_fm_create(work_fm, full_struct)
    1649           8 :       CALL cp_fm_create(sggs_fm, sggs_struct)
    1650             : 
    1651          16 :       ALLOCATE (diag(nactive))
    1652          32 :       ALLOCATE (mo_op(nactive, nactive))
    1653          24 :       ALLOCATE (sggs_block(nactive, nactive))
    1654             : 
    1655             :       ! Iterate over the dimensions of the operator
    1656             :       ! Note: operator matrices are asusmed symmetric, can only do upper half
    1657          32 :       DO i = 1, dim_op
    1658             : 
    1659          24 :          ao_op_i => ao_op(i)%matrix
    1660             : 
    1661             :          ! The GS-GS contribution
    1662          24 :          CALL cp_fm_set_element(amew_op(i), 1, 1, gsgs_op(i))
    1663             : 
    1664             :          ! Compute the operator in the MO basis
    1665          24 :          CALL cp_dbcsr_sm_fm_multiply(ao_op_i, gs_coeffs, vec_op, ncol=nactive)
    1666          24 :          CALL cp_fm_get_submatrix(gs_fm, mo_op)  ! instead of prod_fm
    1667             : 
    1668             :          ! Compute the ground-state/singlet components. ao_op*gs_coeffs already stored in vec_op
    1669          24 :          IF (vel_rep) THEN
    1670             :             CALL dip_vel_op(soc_env, qs_env, evects_sing, dbcsr_prod, i, .FALSE., &
    1671           6 :                             gs_coeffs, sggs_fm)
    1672             :          ELSE
    1673             :             CALL parallel_gemm('T', 'N', nactive*nsg, nactive, nao, &
    1674          18 :                                1.0_dp, soc_env%a_coeff, vec_op, 0.0_dp, sggs_fm)
    1675             :          END IF
    1676             : 
    1677          90 :          DO isg = 1, nsg
    1678             :             CALL cp_fm_get_submatrix(fm=sggs_fm, &
    1679             :                                      target_m=sggs_block, &
    1680             :                                      start_row=(isg - 1)*nactive + 1, &
    1681             :                                      start_col=1, &
    1682             :                                      n_rows=nactive, &
    1683          66 :                                      n_cols=nactive)
    1684          66 :             diag(:) = get_diag(sggs_block)
    1685         528 :             op = sqrt2*SUM(diag)
    1686          90 :             CALL cp_fm_set_element(amew_op(i), 1, 1 + isg, op)
    1687             :          END DO !isg
    1688             : 
    1689             :          ! do the singlet-singlet components
    1690             :          !start with the overlap
    1691             :          CALL dbcsr_multiply('N', 'N', 1.0_dp, &
    1692             :                              matrix_s(1)%matrix, dbcsr_sg, &
    1693             :                              0.0_dp, dbcsr_work, &
    1694          24 :                              filter_eps=eps_filter)
    1695             :          CALL dbcsr_multiply('T', 'N', 1.0_dp, &
    1696             :                              dbcsr_sg, dbcsr_work, &
    1697             :                              0.0_dp, dbcsr_ovlp, &
    1698          24 :                              filter_eps=eps_filter)
    1699             : 
    1700          24 :          IF (vel_rep) THEN
    1701           6 :             CALL dip_vel_op(soc_env, qs_env, evects_sing, dbcsr_prod, i, .FALSE.)
    1702             :          ELSE
    1703             :             !then the operator in the LR orbital basis
    1704             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, &
    1705             :                                 ao_op_i, dbcsr_sg, &
    1706             :                                 0.0_dp, dbcsr_work, &
    1707          18 :                                 filter_eps=eps_filter)
    1708             :             CALL dbcsr_multiply('T', 'N', 1.0_dp, &
    1709             :                                 dbcsr_sg, dbcsr_work, &
    1710             :                                 0.0_dp, dbcsr_prod, &
    1711          18 :                                 filter_eps=eps_filter)
    1712             :          END IF
    1713             : 
    1714             :          !use the soc routine, it is compatible
    1715             :          CALL rcs_amew_soc_elements(dbcsr_tmp, &
    1716             :                                     dbcsr_prod, &
    1717             :                                     dbcsr_ovlp, &
    1718             :                                     mo_op, &
    1719             :                                     pref_trace=-1.0_dp, &
    1720             :                                     pref_overall=1.0_dp, &
    1721             :                                     pref_diags=gsgs_op(i), &
    1722          24 :                                     symmetric=.TRUE.)
    1723             : 
    1724          24 :          CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
    1725             :          CALL cp_fm_to_fm_submat(msource=tmp_fm, &
    1726             :                                  mtarget=amew_op(i), &
    1727             :                                  nrow=nex, &
    1728             :                                  ncol=nex, &
    1729             :                                  s_firstrow=1, &
    1730             :                                  s_firstcol=1, &
    1731             :                                  t_firstrow=2, &
    1732          24 :                                  t_firstcol=2)
    1733             : 
    1734             :          ! compute the triplet-triplet components
    1735             :          !the overlap
    1736             :          CALL dbcsr_multiply('N', 'N', 1.0_dp, &
    1737             :                              matrix_s(1)%matrix, &
    1738             :                              dbcsr_tp, 0.0_dp, &
    1739             :                              dbcsr_work, &
    1740          24 :                              filter_eps=eps_filter)
    1741             :          CALL dbcsr_multiply('T', 'N', 1.0_dp, &
    1742             :                              dbcsr_tp, dbcsr_work, &
    1743             :                              0.0_dp, dbcsr_ovlp, &
    1744          24 :                              filter_eps=eps_filter)
    1745             : 
    1746             :          !the operator in the LR orbital basis
    1747          24 :          IF (vel_rep) THEN
    1748           6 :             CALL dip_vel_op(soc_env, qs_env, evects_trip, dbcsr_prod, i, .TRUE.)
    1749             :          ELSE
    1750             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, &
    1751             :                                 ao_op_i, dbcsr_tp, &
    1752             :                                 0.0_dp, dbcsr_work, &
    1753          18 :                                 filter_eps=eps_filter)
    1754             :             CALL dbcsr_multiply('T', 'N', 1.0_dp, &
    1755             :                                 dbcsr_tp, dbcsr_work, &
    1756             :                                 0.0_dp, dbcsr_prod, &
    1757          18 :                                 filter_eps=eps_filter)
    1758             :          END IF
    1759             : 
    1760             :          CALL rcs_amew_soc_elements(dbcsr_tmp, &
    1761             :                                     dbcsr_prod, &
    1762             :                                     dbcsr_ovlp, &
    1763             :                                     mo_op, &
    1764             :                                     pref_trace=-1.0_dp, &
    1765             :                                     pref_overall=1.0_dp, &
    1766             :                                     pref_diags=gsgs_op(i), &
    1767          24 :                                     symmetric=.TRUE.)
    1768             : 
    1769          24 :          CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
    1770             :          !<T^-1|op|T^-1>
    1771             :          CALL cp_fm_to_fm_submat(msource=tmp_fm, &
    1772             :                                  mtarget=amew_op(i), &
    1773             :                                  nrow=nex, &
    1774             :                                  ncol=nex, &
    1775             :                                  s_firstrow=1, &
    1776             :                                  s_firstcol=1, &
    1777             :                                  t_firstrow=1 + nsg + 1, &
    1778          24 :                                  t_firstcol=1 + nsg + 1)
    1779             : 
    1780             :          !<T^0|op|T^0>
    1781             :          CALL cp_fm_to_fm_submat(msource=tmp_fm, &
    1782             :                                  mtarget=amew_op(i), &
    1783             :                                  nrow=nex, &
    1784             :                                  ncol=nex, &
    1785             :                                  s_firstrow=1, &
    1786             :                                  s_firstcol=1, &
    1787             :                                  t_firstrow=1 + nsg + ntp + 1, &
    1788          24 :                                  t_firstcol=1 + nsg + ntp + 1)
    1789             :          !<T^+1|op|T^+1>
    1790             :          CALL cp_fm_to_fm_submat(msource=tmp_fm, &
    1791             :                                  mtarget=amew_op(i), &
    1792             :                                  nrow=nex, &
    1793             :                                  ncol=nex, &
    1794             :                                  s_firstrow=1, &
    1795             :                                  s_firstcol=1, &
    1796             :                                  t_firstrow=1 + nsg + 2*ntp + 1, &
    1797          24 :                                  t_firstcol=1 + nsg + 2*ntp + 1)
    1798             : 
    1799             :          ! Symmetrize the matrix (only upper triangle built)
    1800          32 :          CALL cp_fm_uplo_to_full(amew_op(i), work_fm)
    1801             : 
    1802             :       END DO !i
    1803             : 
    1804             :       ! Clean-up
    1805           8 :       CALL cp_fm_release(gs_fm)
    1806           8 :       CALL cp_fm_release(work_fm)
    1807           8 :       CALL cp_fm_release(tmp_fm)
    1808           8 :       CALL cp_fm_release(vec_op)
    1809           8 :       CALL cp_fm_release(sggs_fm)
    1810           8 :       CALL cp_fm_struct_release(full_struct)
    1811           8 :       CALL cp_fm_struct_release(tmp_struct)
    1812           8 :       CALL cp_fm_struct_release(sggs_struct)
    1813             : 
    1814           8 :       DEALLOCATE (sggs_block)
    1815           8 :       DEALLOCATE (diag, gsgs_op, mo_op)
    1816           8 :       CALL timestop(handle)
    1817             : 
    1818          40 :    END SUBROUTINE get_rcs_amew_op
    1819             : 
    1820           0 : END MODULE qs_tddfpt2_soc

Generated by: LCOV version 1.15