LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_soc_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 94.3 % 211 199
Test Date: 2025-07-25 12:55:17 Functions: 87.5 % 8 7

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Utilities absorption spectroscopy using TDDFPT with SOC
      10              : !> \author JRVogt (12.2023)
      11              : ! **************************************************************************************************
      12              : 
      13              : MODULE qs_tddfpt2_soc_utils
      14              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      15              :    USE cp_cfm_types,                    ONLY: cp_cfm_get_info,&
      16              :                                               cp_cfm_get_submatrix,&
      17              :                                               cp_cfm_type
      18              :    USE cp_control_types,                ONLY: tddfpt2_control_type
      19              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      20              :                                               dbcsr_create,&
      21              :                                               dbcsr_desymmetrize,&
      22              :                                               dbcsr_get_info,&
      23              :                                               dbcsr_p_type,&
      24              :                                               dbcsr_release,&
      25              :                                               dbcsr_type
      26              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      27              :                                               copy_fm_to_dbcsr,&
      28              :                                               cp_dbcsr_sm_fm_multiply,&
      29              :                                               dbcsr_allocate_matrix_set,&
      30              :                                               dbcsr_deallocate_matrix_set
      31              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_schur_product
      32              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      33              :                                               cp_fm_struct_release,&
      34              :                                               cp_fm_struct_type
      35              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      36              :                                               cp_fm_get_info,&
      37              :                                               cp_fm_release,&
      38              :                                               cp_fm_set_all,&
      39              :                                               cp_fm_to_fm,&
      40              :                                               cp_fm_to_fm_submat,&
      41              :                                               cp_fm_type
      42              :    USE input_constants,                 ONLY: tddfpt_dipole_berry,&
      43              :                                               tddfpt_dipole_length,&
      44              :                                               tddfpt_dipole_velocity
      45              :    USE kinds,                           ONLY: dp
      46              :    USE message_passing,                 ONLY: mp_para_env_type
      47              :    USE moments_utils,                   ONLY: get_reference_point
      48              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      49              :    USE qs_environment_types,            ONLY: get_qs_env,&
      50              :                                               qs_environment_type
      51              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      52              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      53              :    USE qs_operators_ao,                 ONLY: p_xyz_ao,&
      54              :                                               rRc_xyz_ao
      55              :    USE qs_overlap,                      ONLY: build_overlap_matrix
      56              :    USE qs_tddfpt2_soc_types,            ONLY: soc_env_type
      57              :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos
      58              : 
      59              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
      60              : #include "./base/base_uses.f90"
      61              : 
      62              :    IMPLICIT NONE
      63              :    PRIVATE
      64              : 
      65              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_soc_utils'
      66              : 
      67              :    PUBLIC :: soc_dipole_operator, soc_contract_evect, resort_evects, dip_vel_op
      68              : 
      69              :    !A helper type for SOC
      70              :    TYPE dbcsr_soc_package_type
      71              :       TYPE(dbcsr_type), POINTER     :: dbcsr_sg => Null()
      72              :       TYPE(dbcsr_type), POINTER     :: dbcsr_tp => Null()
      73              :       TYPE(dbcsr_type), POINTER     :: dbcsr_sc => Null()
      74              :       TYPE(dbcsr_type), POINTER     :: dbcsr_sf => Null()
      75              :       TYPE(dbcsr_type), POINTER     :: dbcsr_prod => Null()
      76              :       TYPE(dbcsr_type), POINTER     :: dbcsr_ovlp => Null()
      77              :       TYPE(dbcsr_type), POINTER     :: dbcsr_tmp => Null()
      78              :       TYPE(dbcsr_type), POINTER     :: dbcsr_work => Null()
      79              :    END TYPE dbcsr_soc_package_type
      80              : 
      81              : CONTAINS
      82              : 
      83              : ! **************************************************************************************************
      84              : !> \brief Build the atomic dipole operator
      85              : !> \param soc_env ...
      86              : !> \param tddfpt_control informations on how to build the operaot
      87              : !> \param qs_env Qucikstep environment
      88              : !> \param gs_mos ...
      89              : ! **************************************************************************************************
      90            8 :    SUBROUTINE soc_dipole_operator(soc_env, tddfpt_control, qs_env, gs_mos)
      91              :       TYPE(soc_env_type), TARGET                         :: soc_env
      92              :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
      93              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
      94              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
      95              :          INTENT(in)                                      :: gs_mos
      96              : 
      97              :       CHARACTER(len=*), PARAMETER :: routineN = 'soc_dipole_operator'
      98              : 
      99              :       INTEGER                                            :: dim_op, handle, i_dim, nao, nspin
     100              :       REAL(kind=dp), DIMENSION(3)                        :: reference_point
     101            8 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     102              : 
     103            8 :       CALL timeset(routineN, handle)
     104              : 
     105            8 :       NULLIFY (matrix_s)
     106              : 
     107            8 :       IF (tddfpt_control%dipole_form == tddfpt_dipole_berry) THEN
     108            0 :          CPABORT("BERRY DIPOLE FORM NOT IMPLEMENTED FOR SOC")
     109              :       END IF
     110              :            !! ONLY RCS have been implemented, Therefore, nspin sould always be 1!
     111            8 :       nspin = 1
     112              :            !! Number of dimensions should be 3, unless multipole is implemented in the future
     113            8 :       dim_op = 3
     114              : 
     115              :            !! Initzilize the dipmat structure
     116            8 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     117            8 :       CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
     118              : 
     119           32 :       ALLOCATE (soc_env%dipmat_ao(dim_op))
     120           32 :       DO i_dim = 1, dim_op
     121           24 :          ALLOCATE (soc_env%dipmat_ao(i_dim)%matrix)
     122              :          CALL dbcsr_copy(soc_env%dipmat_ao(i_dim)%matrix, &
     123              :                          matrix_s(1)%matrix, &
     124           32 :                          name="dipole operator matrix")
     125              :       END DO
     126              : 
     127           14 :       SELECT CASE (tddfpt_control%dipole_form)
     128              :       CASE (tddfpt_dipole_length)
     129              :                    !!This routine is analog to qs_tddfpt_prperties but only until the rRc_xyz_ao routine
     130              :                    !! This will lead to an operator within the nao x nao basis
     131              :                    !! qs_tddpft_properies uses nvirt x nocc
     132              :          CALL get_reference_point(reference_point, qs_env=qs_env, &
     133              :                                   reference=tddfpt_control%dipole_reference, &
     134            6 :                                   ref_point=tddfpt_control%dipole_ref_point)
     135              : 
     136              :          CALL rRc_xyz_ao(op=soc_env%dipmat_ao, qs_env=qs_env, rc=reference_point, order=1, &
     137            6 :                          minimum_image=.FALSE., soft=.FALSE.)
     138              :          !! This will lead to S C^virt C^virt,T Q_q (vgl Strand et al., J. Chem Phys. 150, 044702, 2019)
     139            6 :          CALL length_rep(qs_env, gs_mos, soc_env)
     140              :       CASE (tddfpt_dipole_velocity)
     141              :          !!This Routine calcluates the dipole Operator within the velocity-form within the ao basis
     142              :          !!This Operation is only used in xas_tdp and qs_tddfpt_soc, lines uses rmc_x_p_xyz_ao
     143            2 :          CALL p_xyz_ao(soc_env%dipmat_ao, qs_env, minimum_image=.FALSE.)
     144              :          !! This will precomute SC^virt, (omega^a-omega^i)^-1 and C^virt dS/dq
     145            2 :          CALL velocity_rep(qs_env, gs_mos, soc_env)
     146              :       CASE DEFAULT
     147            8 :          CPABORT("Unimplemented form of the dipole operator")
     148              :       END SELECT
     149              : 
     150            8 :       CALL timestop(handle)
     151              : 
     152           16 :    END SUBROUTINE soc_dipole_operator
     153              : 
     154              : ! **************************************************************************************************
     155              : !> \brief ...
     156              : !> \param qs_env ...
     157              : !> \param gs_mos ...
     158              : !> \param soc_env ...
     159              : ! **************************************************************************************************
     160            6 :    SUBROUTINE length_rep(qs_env, gs_mos, soc_env)
     161              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     162              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     163              :          INTENT(in)                                      :: gs_mos
     164              :       TYPE(soc_env_type), TARGET                         :: soc_env
     165              : 
     166              :       INTEGER                                            :: ideriv, ispin, nao, nderivs, nspins
     167            6 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nmo_virt
     168              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     169              :       TYPE(cp_fm_struct_type), POINTER                   :: dip_struct, fm_struct
     170            6 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: S_mos_virt
     171            6 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: dipole_op_mos_occ
     172              :       TYPE(cp_fm_type), POINTER                          :: dipmat_tmp, wfm_ao_ao
     173            6 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     174              :       TYPE(dbcsr_type), POINTER                          :: symm_tmp
     175              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     176              : 
     177            6 :       CALL get_qs_env(qs_env, matrix_s=matrix_s, blacs_env=blacs_env, para_env=para_env)
     178              : 
     179            6 :       nderivs = 3
     180            6 :       nspins = 1  !!We only account for rcs, will be changed in the future
     181            6 :       CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
     182              :       ALLOCATE (S_mos_virt(nspins), dipole_op_mos_occ(3, nspins), &
     183           36 :                 wfm_ao_ao, nmo_virt(nspins), symm_tmp, dipmat_tmp)
     184              : 
     185            6 :       CALL cp_fm_struct_create(dip_struct, context=blacs_env, ncol_global=nao, nrow_global=nao, para_env=para_env)
     186              : 
     187            6 :       CALL dbcsr_allocate_matrix_set(soc_env%dipmat, nderivs)
     188            6 :       CALL dbcsr_desymmetrize(matrix_s(1)%matrix, symm_tmp)
     189           24 :       DO ideriv = 1, nderivs
     190           18 :          ALLOCATE (soc_env%dipmat(ideriv)%matrix)
     191              :          CALL dbcsr_create(soc_env%dipmat(ideriv)%matrix, template=symm_tmp, &
     192           18 :                            name="contracted operator", matrix_type="N")
     193           42 :          DO ispin = 1, nspins
     194           36 :             CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin), matrix_struct=dip_struct)
     195              :          END DO
     196              :       END DO
     197              : 
     198            6 :       CALL dbcsr_release(symm_tmp)
     199            6 :       DEALLOCATE (symm_tmp)
     200              : 
     201           12 :       DO ispin = 1, nspins
     202            6 :          nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
     203            6 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
     204            6 :          CALL cp_fm_create(wfm_ao_ao, dip_struct)
     205            6 :          CALL cp_fm_create(S_mos_virt(ispin), fm_struct)
     206              : 
     207              :          CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
     208              :                                       gs_mos(ispin)%mos_virt, &
     209              :                                       S_mos_virt(ispin), &
     210            6 :                                       ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
     211              :          CALL parallel_gemm('N', 'T', nao, nao, nmo_virt(ispin), &
     212              :                             1.0_dp, S_mos_virt(ispin), gs_mos(ispin)%mos_virt, &
     213            6 :                             0.0_dp, wfm_ao_ao)
     214              : 
     215           24 :          DO ideriv = 1, nderivs
     216           18 :             CALL cp_fm_create(dipmat_tmp, dip_struct)
     217           18 :             CALL copy_dbcsr_to_fm(soc_env%dipmat_ao(ideriv)%matrix, dipmat_tmp)
     218              :             CALL parallel_gemm('N', 'T', nao, nao, nao, &
     219              :                                1.0_dp, wfm_ao_ao, dipmat_tmp, &
     220           18 :                                0.0_dp, dipole_op_mos_occ(ideriv, ispin))
     221           18 :             CALL copy_fm_to_dbcsr(dipole_op_mos_occ(ideriv, ispin), soc_env%dipmat(ideriv)%matrix)
     222           42 :             CALL cp_fm_release(dipmat_tmp)
     223              :          END DO
     224            6 :          CALL cp_fm_release(wfm_ao_ao)
     225           18 :          DEALLOCATE (wfm_ao_ao)
     226              :       END DO
     227              : 
     228            6 :       CALL cp_fm_struct_release(dip_struct)
     229           12 :       DO ispin = 1, nspins
     230            6 :          CALL cp_fm_release(S_mos_virt(ispin))
     231           30 :          DO ideriv = 1, nderivs
     232           24 :             CALL cp_fm_release(dipole_op_mos_occ(ideriv, ispin))
     233              :          END DO
     234              :       END DO
     235            6 :       DEALLOCATE (S_mos_virt, dipole_op_mos_occ, nmo_virt, dipmat_tmp)
     236              : 
     237            6 :    END SUBROUTINE length_rep
     238              : 
     239              : ! **************************************************************************************************
     240              : !> \brief ...
     241              : !> \param qs_env ...
     242              : !> \param gs_mos ...
     243              : !> \param soc_env ...
     244              : ! **************************************************************************************************
     245            2 :    SUBROUTINE velocity_rep(qs_env, gs_mos, soc_env)
     246              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     247              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     248              :          INTENT(in)                                      :: gs_mos
     249              :       TYPE(soc_env_type), TARGET                         :: soc_env
     250              : 
     251              :       INTEGER                                            :: icol, ideriv, irow, ispin, n_occ, &
     252              :                                                             n_virt, nao, ncols_local, nderivs, &
     253              :                                                             nrows_local, nspins
     254            2 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     255              :       REAL(kind=dp)                                      :: eval_occ
     256              :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
     257            2 :          POINTER                                         :: local_data_ediff, local_data_wfm
     258              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     259              :       TYPE(cp_fm_struct_type), POINTER                   :: ao_cvirt_struct, ao_nocc_struct, &
     260              :                                                             cvirt_ao_struct, fm_struct, scrm_struct
     261              :       TYPE(cp_fm_type)                                   :: scrm_fm, wfm_mo_virt_mo_occ
     262            2 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, scrm
     263              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     264            2 :          POINTER                                         :: sab_orb
     265              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     266              : 
     267            2 :       NULLIFY (scrm, scrm_struct, blacs_env, matrix_s, ao_cvirt_struct, ao_nocc_struct, cvirt_ao_struct)
     268            2 :       nspins = 1
     269            2 :       nderivs = 3
     270           18 :       ALLOCATE (soc_env%SC(nspins), soc_env%CdS(nspins, nderivs), soc_env%ediff(nspins))
     271              : 
     272            2 :       CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb, blacs_env=blacs_env, matrix_s=matrix_s)
     273            2 :       CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
     274              :       CALL cp_fm_struct_create(scrm_struct, nrow_global=nao, ncol_global=nao, &
     275            2 :                                context=blacs_env)
     276            2 :       CALL cp_fm_get_info(gs_mos(1)%mos_virt, matrix_struct=ao_cvirt_struct)
     277            2 :       CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=ao_nocc_struct)
     278              : 
     279              :       CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
     280              :                                 basis_type_a="ORB", basis_type_b="ORB", &
     281            2 :                                 sab_nl=sab_orb)
     282              : 
     283            4 :       DO ispin = 1, nspins
     284            2 :          NULLIFY (fm_struct)
     285            2 :          n_occ = SIZE(gs_mos(ispin)%evals_occ)
     286            2 :          n_virt = SIZE(gs_mos(ispin)%evals_virt)
     287              :          CALL cp_fm_struct_create(fm_struct, nrow_global=n_virt, &
     288            2 :                                   ncol_global=n_occ, context=blacs_env)
     289              :          CALL cp_fm_struct_create(cvirt_ao_struct, nrow_global=n_virt, &
     290            2 :                                   ncol_global=nao, context=blacs_env)
     291            2 :          CALL cp_fm_create(soc_env%ediff(ispin), fm_struct)
     292            2 :          CALL cp_fm_create(wfm_mo_virt_mo_occ, fm_struct)
     293            2 :          CALL cp_fm_create(soc_env%SC(ispin), ao_cvirt_struct)
     294              : 
     295              :          CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
     296              :                                       gs_mos(ispin)%mos_virt, &
     297              :                                       soc_env%SC(ispin), &
     298            2 :                                       ncol=n_virt, alpha=1.0_dp, beta=0.0_dp)
     299              : 
     300              :          CALL cp_fm_get_info(soc_env%ediff(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
     301            2 :                              row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
     302            2 :          CALL cp_fm_get_info(wfm_mo_virt_mo_occ, local_data=local_data_wfm)
     303              : 
     304              : !$OMP       PARALLEL DO DEFAULT(NONE), &
     305              : !$OMP                PRIVATE(eval_occ, icol, irow), &
     306            2 : !$OMP                SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices)
     307              :          DO icol = 1, ncols_local
     308              :             ! E_occ_i ; imo_occ = col_indices(icol)
     309              :             eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol))
     310              : 
     311              :             DO irow = 1, nrows_local
     312              :                ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i)
     313              :                ! imo_virt = row_indices(irow)
     314              :                local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ)
     315              :             END DO
     316              :          END DO
     317              : !$OMP       END PARALLEL DO
     318              : 
     319            8 :          DO ideriv = 1, nderivs
     320            6 :             CALL cp_fm_create(soc_env%CdS(ispin, ideriv), cvirt_ao_struct)
     321            6 :             CALL cp_fm_create(scrm_fm, scrm_struct)
     322            6 :             CALL copy_dbcsr_to_fm(scrm(ideriv + 1)%matrix, scrm_fm)
     323              :             CALL parallel_gemm('T', 'N', n_virt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
     324            6 :                                scrm_fm, 0.0_dp, soc_env%CdS(ispin, ideriv))
     325           14 :             CALL cp_fm_release(scrm_fm)
     326              : 
     327              :          END DO
     328              : 
     329            2 :          CALL cp_fm_release(wfm_mo_virt_mo_occ)
     330            8 :          CALL cp_fm_struct_release(fm_struct)
     331              :       END DO
     332            2 :       CALL dbcsr_deallocate_matrix_set(scrm)
     333            2 :       CALL cp_fm_struct_release(scrm_struct)
     334            2 :       CALL cp_fm_struct_release(cvirt_ao_struct)
     335              : 
     336            4 :    END SUBROUTINE velocity_rep
     337              : 
     338              : ! **************************************************************************************************
     339              : !> \brief This routine will construct the dipol operator within velocity representation
     340              : !> \param soc_env ..
     341              : !> \param qs_env ...
     342              : !> \param evec_fm ...
     343              : !> \param op ...
     344              : !> \param ideriv ...
     345              : !> \param tp ...
     346              : !> \param gs_coeffs ...
     347              : !> \param sggs_fm ...
     348              : ! **************************************************************************************************
     349           18 :    SUBROUTINE dip_vel_op(soc_env, qs_env, evec_fm, op, ideriv, tp, gs_coeffs, sggs_fm)
     350              :       TYPE(soc_env_type), TARGET                         :: soc_env
     351              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     352              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: evec_fm
     353              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: op
     354              :       INTEGER, INTENT(IN)                                :: ideriv
     355              :       LOGICAL, INTENT(IN)                                :: tp
     356              :       TYPE(cp_fm_type), OPTIONAL, POINTER                :: gs_coeffs
     357              :       TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL          :: sggs_fm
     358              : 
     359              :       INTEGER                                            :: iex, ispin, n_occ, n_virt, nao, nex
     360              :       LOGICAL                                            :: sggs
     361              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     362              :       TYPE(cp_fm_struct_type), POINTER                   :: op_struct, virt_occ_struct
     363              :       TYPE(cp_fm_type)                                   :: CdSC, op_fm, SCWCdSC, WCdSC
     364           18 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: WCdSC_tmp
     365              :       TYPE(cp_fm_type), POINTER                          :: coeff
     366              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     367              : 
     368           18 :       NULLIFY (virt_occ_struct, virt_occ_struct, op_struct, blacs_env, para_env, coeff)
     369              : 
     370           18 :       IF (tp) THEN
     371            6 :          coeff => soc_env%b_coeff
     372              :       ELSE
     373           12 :          coeff => soc_env%a_coeff
     374              :       END IF
     375              : 
     376           18 :       sggs = .FALSE.
     377           18 :       IF (PRESENT(gs_coeffs)) sggs = .TRUE.
     378              : 
     379           18 :       ispin = 1 !! only rcs availble
     380           18 :       nex = SIZE(evec_fm, 2)
     381           90 :       IF (.NOT. sggs) ALLOCATE (WCdSC_tmp(ispin, nex))
     382           18 :       CALL get_qs_env(qs_env, blacs_env=blacs_env, para_env=para_env)
     383           18 :       CALL cp_fm_get_info(soc_env%CdS(ispin, ideriv), ncol_global=nao, nrow_global=n_virt)
     384           18 :       CALL cp_fm_get_info(evec_fm(1, 1), ncol_global=n_occ)
     385              : 
     386           18 :       IF (sggs) THEN
     387              :          CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
     388            6 :                                   ncol_global=n_occ)
     389              :          CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_occ*nex, &
     390            6 :                                   ncol_global=n_occ)
     391              :       ELSE
     392              :          CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
     393           12 :                                   ncol_global=n_occ*nex)
     394              :          CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_occ*nex, &
     395           12 :                                   ncol_global=n_occ*nex)
     396              :       END IF
     397              : 
     398           18 :       CALL cp_fm_create(CdSC, soc_env%ediff(ispin)%matrix_struct)
     399           18 :       CALL cp_fm_create(op_fm, op_struct)
     400              : 
     401           18 :       IF (sggs) THEN
     402            6 :          CALL cp_fm_create(SCWCdSC, gs_coeffs%matrix_struct)
     403            6 :          CALL cp_fm_create(WCdSC, soc_env%ediff(ispin)%matrix_struct)
     404              :          CALL parallel_gemm('N', 'N', n_virt, n_occ, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
     405            6 :                             gs_coeffs, 0.0_dp, CdSC)
     406            6 :          CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC)
     407              :       ELSE
     408           12 :          CALL cp_fm_create(SCWCdSC, coeff%matrix_struct)
     409           36 :          DO iex = 1, nex
     410           24 :             CALL cp_fm_create(WCdSC_tmp(ispin, iex), soc_env%ediff(ispin)%matrix_struct)
     411              :             CALL parallel_gemm('N', 'N', n_virt, n_occ, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
     412           24 :                                evec_fm(ispin, iex), 0.0_dp, CdSC)
     413           36 :             CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC_tmp(ispin, iex))
     414              :          END DO
     415           12 :          CALL cp_fm_create(WCdSC, virt_occ_struct)
     416           12 :          CALL soc_contract_evect(WCdSC_tmp, WCdSC)
     417           36 :          DO iex = 1, nex
     418           36 :             CALL cp_fm_release(WCdSC_tmp(ispin, iex))
     419              :          END DO
     420           12 :          DEALLOCATE (WCdSC_tmp)
     421              :       END IF
     422              : 
     423           18 :       IF (sggs) THEN
     424            6 :          CALL parallel_gemm('N', 'N', nao, n_occ, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
     425            6 :          CALL parallel_gemm('T', 'N', n_occ*nex, n_occ, nao, 1.0_dp, soc_env%a_coeff, SCWCdSC, 0.0_dp, op_fm)
     426              :       ELSE
     427           12 :          CALL parallel_gemm('N', 'N', nao, n_occ*nex, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
     428           12 :          CALL parallel_gemm('T', 'N', n_occ*nex, n_occ*nex, nao, 1.0_dp, coeff, SCWCdSC, 0.0_dp, op_fm)
     429              :       END IF
     430              : 
     431           18 :       IF (sggs) THEN
     432            6 :          CALL cp_fm_to_fm(op_fm, sggs_fm)
     433              :       ELSE
     434           12 :          CALL copy_fm_to_dbcsr(op_fm, op)
     435              :       END IF
     436              : 
     437           18 :       CALL cp_fm_release(op_fm)
     438           18 :       CALL cp_fm_release(WCdSC)
     439           18 :       CALL cp_fm_release(SCWCdSC)
     440           18 :       CALL cp_fm_release(CdSC)
     441           18 :       CALL cp_fm_struct_release(virt_occ_struct)
     442           18 :       CALL cp_fm_struct_release(op_struct)
     443              : 
     444           36 :    END SUBROUTINE dip_vel_op
     445              : 
     446              : ! **************************************************************************************************
     447              : !> \brief ...
     448              : !> \param fm_start ...
     449              : !> \param fm_res ...
     450              : ! **************************************************************************************************
     451           28 :    SUBROUTINE soc_contract_evect(fm_start, fm_res)
     452              : 
     453              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: fm_start
     454              :       TYPE(cp_fm_type), INTENT(inout)                    :: fm_res
     455              : 
     456              :       CHARACTER(len=*), PARAMETER :: routineN = 'soc_contract_evect'
     457              : 
     458              :       INTEGER                                            :: handle, ii, jj, nactive, nao, nspins, &
     459              :                                                             nstates, ntmp1, ntmp2
     460              : 
     461           28 :       CALL timeset(routineN, handle)
     462              : 
     463           28 :       nstates = SIZE(fm_start, 2)
     464           28 :       nspins = SIZE(fm_start, 1)
     465              : 
     466           28 :       CALL cp_fm_set_all(fm_res, 0.0_dp)
     467              :          !! Evects are written into one matrix.
     468           96 :       DO ii = 1, nstates
     469          164 :          DO jj = 1, nspins
     470           68 :             CALL cp_fm_get_info(fm_start(jj, ii), nrow_global=nao, ncol_global=nactive)
     471           68 :             CALL cp_fm_get_info(fm_res, nrow_global=ntmp1, ncol_global=ntmp2)
     472              :             CALL cp_fm_to_fm_submat(fm_start(jj, ii), &
     473              :                                     fm_res, &
     474              :                                     nao, nactive, &
     475              :                                     1, 1, 1, &
     476          136 :                                     1 + nactive*(ii - 1) + (jj - 1)*nao*nstates)
     477              :          END DO !nspins
     478              :       END DO !nsstates
     479              : 
     480           28 :       CALL timestop(handle)
     481              : 
     482           28 :    END SUBROUTINE soc_contract_evect
     483              : 
     484              : ! **************************************************************************************************
     485              : !> \brief ...
     486              : !> \param vec ...
     487              : !> \param new_entry ...
     488              : !> \param res ...
     489              : !> \param res_int ...
     490              : ! **************************************************************************************************
     491          454 :    SUBROUTINE test_repetition(vec, new_entry, res, res_int)
     492              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: vec
     493              :       INTEGER, INTENT(IN)                                :: new_entry
     494              :       LOGICAL, INTENT(OUT)                               :: res
     495              :       INTEGER, INTENT(OUT), OPTIONAL                     :: res_int
     496              : 
     497              :       INTEGER                                            :: i
     498              : 
     499          454 :       res = .TRUE.
     500          454 :       IF (PRESENT(res_int)) res_int = -1
     501              : 
     502         4432 :       DO i = 1, SIZE(vec)
     503         4432 :          IF (vec(i) == new_entry) THEN
     504          166 :             res = .FALSE.
     505          166 :             IF (PRESENT(res_int)) res_int = i
     506              :             EXIT
     507              :          END IF
     508              :       END DO
     509              : 
     510          454 :    END SUBROUTINE test_repetition
     511              : 
     512              : ! **************************************************************************************************
     513              : !> \brief Used to find out, which state has which spin-multiplicity
     514              : !> \param evects_cfm ...
     515              : !> \param sort ...
     516              : ! **************************************************************************************************
     517            8 :    SUBROUTINE resort_evects(evects_cfm, sort)
     518              :       TYPE(cp_cfm_type), INTENT(INOUT)                   :: evects_cfm
     519              :       INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT)    :: sort
     520              : 
     521              :       COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :)          :: cpl_tmp
     522              :       INTEGER                                            :: i_rep, ii, jj, ntot, tmp
     523            8 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: rep_int
     524              :       LOGICAL                                            :: rep
     525              :       REAL(dp)                                           :: max_dev, max_wfn, wfn_sq
     526              : 
     527            8 :       CALL cp_cfm_get_info(evects_cfm, nrow_global=ntot)
     528           32 :       ALLOCATE (cpl_tmp(ntot, ntot))
     529           32 :       ALLOCATE (sort(ntot), rep_int(ntot))
     530         1280 :       cpl_tmp = 0_dp
     531          104 :       sort = 0
     532            8 :       max_dev = 0.5
     533            8 :       CALL cp_cfm_get_submatrix(evects_cfm, cpl_tmp)
     534              : 
     535          104 :       DO jj = 1, ntot
     536         1272 :          rep_int = 0
     537           96 :          tmp = 0
     538           96 :          max_wfn = 0_dp
     539         1272 :          DO ii = 1, ntot
     540         1176 :             wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
     541         1272 :             IF (max_wfn .LE. wfn_sq) THEN
     542          454 :                CALL test_repetition(sort, ii, rep, rep_int(ii))
     543          454 :                IF (rep) THEN
     544          288 :                   max_wfn = wfn_sq
     545          288 :                   tmp = ii
     546              :                END IF
     547              :             END IF
     548              :          END DO
     549          104 :          IF (tmp > 0) THEN
     550           96 :             sort(jj) = tmp
     551              :          ELSE
     552            0 :             DO i_rep = 1, ntot
     553            0 :                IF (rep_int(i_rep) > 0) THEN
     554            0 :                   max_wfn = ABS(REAL(cpl_tmp(sort(i_rep), jj)**2 - AIMAG(cpl_tmp(sort(i_rep), jj)**2))) - max_dev
     555            0 :                   DO ii = 1, ntot
     556            0 :                      wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
     557            0 :                      IF ((max_wfn - wfn_sq)/max_wfn .LE. max_dev) THEN
     558            0 :                         CALL test_repetition(sort, ii, rep)
     559            0 :                         IF (rep .AND. ii /= i_rep) THEN
     560            0 :                            sort(jj) = sort(i_rep)
     561            0 :                            sort(i_rep) = ii
     562              :                         END IF
     563              :                      END IF
     564              :                   END DO
     565              :                END IF
     566              :             END DO
     567              :          END IF
     568              :       END DO
     569              : 
     570            8 :       DEALLOCATE (cpl_tmp, rep_int)
     571              : 
     572            8 :    END SUBROUTINE resort_evects
     573            0 : END MODULE qs_tddfpt2_soc_utils
        

Generated by: LCOV version 2.0-1