LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_properties.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.1 % 584 538
Test Date: 2025-07-25 12:55:17 Functions: 75.0 % 8 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : MODULE qs_tddfpt2_properties
       9              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      10              :    USE bibliography,                    ONLY: Martin2003,&
      11              :                                               cite_reference
      12              :    USE bse_print,                       ONLY: print_exciton_descriptors
      13              :    USE bse_properties,                  ONLY: exciton_descr_type,&
      14              :                                               get_exciton_descriptors
      15              :    USE bse_util,                        ONLY: get_multipoles_mo
      16              :    USE cell_types,                      ONLY: cell_type
      17              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      18              :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_solve
      19              :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      20              :                                               cp_cfm_release,&
      21              :                                               cp_cfm_set_all,&
      22              :                                               cp_cfm_to_fm,&
      23              :                                               cp_cfm_type,&
      24              :                                               cp_fm_to_cfm
      25              :    USE cp_control_types,                ONLY: dft_control_type,&
      26              :                                               tddfpt2_control_type
      27              :    USE cp_dbcsr_api,                    ONLY: &
      28              :         dbcsr_copy, dbcsr_get_block_p, dbcsr_get_info, dbcsr_init_p, dbcsr_iterator_blocks_left, &
      29              :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      30              :         dbcsr_p_type, dbcsr_set, dbcsr_type
      31              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      32              :                                               copy_fm_to_dbcsr,&
      33              :                                               cp_dbcsr_sm_fm_multiply,&
      34              :                                               dbcsr_allocate_matrix_set,&
      35              :                                               dbcsr_deallocate_matrix_set
      36              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale,&
      37              :                                               cp_fm_scale_and_add,&
      38              :                                               cp_fm_trace
      39              :    USE cp_fm_diag,                      ONLY: choose_eigv_solver,&
      40              :                                               cp_fm_geeig
      41              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      42              :                                               cp_fm_struct_release,&
      43              :                                               cp_fm_struct_type
      44              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      45              :                                               cp_fm_get_info,&
      46              :                                               cp_fm_release,&
      47              :                                               cp_fm_set_all,&
      48              :                                               cp_fm_to_fm,&
      49              :                                               cp_fm_to_fm_submat_general,&
      50              :                                               cp_fm_type,&
      51              :                                               cp_fm_vectorsnorm
      52              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      53              :                                               cp_logger_get_default_io_unit,&
      54              :                                               cp_logger_type
      55              :    USE cp_output_handling,              ONLY: cp_p_file,&
      56              :                                               cp_print_key_finished_output,&
      57              :                                               cp_print_key_should_output,&
      58              :                                               cp_print_key_unit_nr
      59              :    USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
      60              :    USE input_constants,                 ONLY: tddfpt_dipole_berry,&
      61              :                                               tddfpt_dipole_length,&
      62              :                                               tddfpt_dipole_velocity
      63              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      64              :                                               section_vals_type,&
      65              :                                               section_vals_val_get
      66              :    USE kahan_sum,                       ONLY: accurate_dot_product
      67              :    USE kinds,                           ONLY: default_path_length,&
      68              :                                               dp,&
      69              :                                               int_8
      70              :    USE mathconstants,                   ONLY: twopi,&
      71              :                                               z_one,&
      72              :                                               z_zero
      73              :    USE message_passing,                 ONLY: mp_comm_type,&
      74              :                                               mp_para_env_type,&
      75              :                                               mp_request_type
      76              :    USE molden_utils,                    ONLY: write_mos_molden
      77              :    USE moments_utils,                   ONLY: get_reference_point
      78              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      79              :    USE particle_list_types,             ONLY: particle_list_type
      80              :    USE particle_types,                  ONLY: particle_type
      81              :    USE physcon,                         ONLY: evolt
      82              :    USE pw_env_types,                    ONLY: pw_env_get,&
      83              :                                               pw_env_type
      84              :    USE pw_poisson_types,                ONLY: pw_poisson_type
      85              :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      86              :                                               pw_pool_type
      87              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      88              :                                               pw_r3d_rs_type
      89              :    USE qs_collocate_density,            ONLY: calculate_wavefunction
      90              :    USE qs_environment_types,            ONLY: get_qs_env,&
      91              :                                               qs_environment_type
      92              :    USE qs_kind_types,                   ONLY: qs_kind_type
      93              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      94              :    USE qs_mo_types,                     ONLY: allocate_mo_set,&
      95              :                                               deallocate_mo_set,&
      96              :                                               get_mo_set,&
      97              :                                               init_mo_set,&
      98              :                                               mo_set_type,&
      99              :                                               set_mo_set
     100              :    USE qs_moments,                      ONLY: build_berry_moment_matrix
     101              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
     102              :    USE qs_operators_ao,                 ONLY: rRc_xyz_ao
     103              :    USE qs_overlap,                      ONLY: build_overlap_matrix
     104              :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
     105              :                                               qs_subsys_type
     106              :    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos
     107              :    USE string_utilities,                ONLY: integer_to_string
     108              :    USE util,                            ONLY: sort
     109              : #include "./base/base_uses.f90"
     110              : 
     111              :    IMPLICIT NONE
     112              : 
     113              :    PRIVATE
     114              : 
     115              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_properties'
     116              : 
     117              :    ! number of first derivative components (3: d/dx, d/dy, d/dz)
     118              :    INTEGER, PARAMETER, PRIVATE          :: nderivs = 3
     119              :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
     120              : 
     121              :    PUBLIC :: tddfpt_dipole_operator, tddfpt_print_summary, tddfpt_print_excitation_analysis, &
     122              :              tddfpt_print_nto_analysis, tddfpt_print_exciton_descriptors
     123              : 
     124              : ! **************************************************************************************************
     125              : 
     126              : CONTAINS
     127              : 
     128              : ! **************************************************************************************************
     129              : !> \brief Compute the action of the dipole operator on the ground state wave function.
     130              : !> \param dipole_op_mos_occ  2-D array [x,y,z ; spin] of matrices where to put the computed quantity
     131              : !>                           (allocated and initialised on exit)
     132              : !> \param tddfpt_control     TDDFPT control parameters
     133              : !> \param gs_mos             molecular orbitals optimised for the ground state
     134              : !> \param qs_env             Quickstep environment
     135              : !> \par History
     136              : !>    * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
     137              : !>    * 06.2018 dipole operator based on the Berry-phase formula [Sergey Chulkov]
     138              : !>    * 08.2018 splited of from 'tddfpt_print_summary' and merged with code from 'tddfpt'
     139              : !>              [Sergey Chulkov]
     140              : !> \note \parblock
     141              : !>       Adapted version of the subroutine find_contributions() which was originally created
     142              : !>       by Thomas Chassaing on 02.2005.
     143              : !>
     144              : !>       The relation between dipole integrals in velocity and length forms are the following:
     145              : !>       \f[<\psi_i|\nabla|\psi_a> = <\psi_i|\vec{r}|\hat{H}\psi_a> - <\hat{H}\psi_i|\vec{r}|\psi_a>
     146              : !>                                 = (\epsilon_a - \epsilon_i) <\psi_i|\vec{r}|\psi_a> .\f],
     147              : !>       due to the commutation identity:
     148              : !>       \f[\vec{r}\hat{H} - \hat{H}\vec{r} = [\vec{r},\hat{H}] = [\vec{r},-1/2 \nabla^2] = \nabla\f] .
     149              : !>       \endparblock
     150              : ! **************************************************************************************************
     151         1076 :    SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_env)
     152              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :), &
     153              :          INTENT(inout)                                   :: dipole_op_mos_occ
     154              :       TYPE(tddfpt2_control_type), POINTER                :: tddfpt_control
     155              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     156              :          INTENT(in)                                      :: gs_mos
     157              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     158              : 
     159              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_dipole_operator'
     160              : 
     161              :       INTEGER                                            :: handle, i_cos_sin, icol, ideriv, irow, &
     162              :                                                             ispin, jderiv, nao, ncols_local, &
     163              :                                                             ndim_periodic, nrows_local, nspins
     164         1076 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     165              :       INTEGER, DIMENSION(maxspins)                       :: nmo_occ, nmo_virt
     166              :       REAL(kind=dp)                                      :: eval_occ
     167              :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
     168         1076 :          POINTER                                         :: local_data_ediff, local_data_wfm
     169              :       REAL(kind=dp), DIMENSION(3)                        :: kvec, reference_point
     170              :       TYPE(cell_type), POINTER                           :: cell
     171              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     172         1076 :       TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: gamma_00, gamma_inv_00
     173              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     174              :       TYPE(cp_fm_type)                                   :: ediff_inv, rRc_mos_occ, wfm_ao_ao, &
     175              :                                                             wfm_mo_virt_mo_occ
     176         1076 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: S_mos_virt
     177         1076 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: dBerry_mos_occ, gamma_real_imag, opvec
     178         1076 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: berry_cossin_xyz, matrix_s, rRc_xyz, scrm
     179              :       TYPE(dft_control_type), POINTER                    :: dft_control
     180              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     181         1076 :          POINTER                                         :: sab_orb
     182              :       TYPE(pw_env_type), POINTER                         :: pw_env
     183              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     184              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     185              : 
     186         1076 :       CALL timeset(routineN, handle)
     187              : 
     188         1076 :       NULLIFY (blacs_env, cell, matrix_s, pw_env)
     189         1076 :       CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, matrix_s=matrix_s, pw_env=pw_env)
     190              : 
     191         1076 :       nspins = SIZE(gs_mos)
     192         1076 :       CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
     193         2276 :       DO ispin = 1, nspins
     194         1200 :          nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
     195         2276 :          nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
     196              :       END DO
     197              : 
     198              :       ! +++ allocate dipole operator matrices (must be deallocated elsewhere)
     199         8028 :       ALLOCATE (dipole_op_mos_occ(nderivs, nspins))
     200         2276 :       DO ispin = 1, nspins
     201         1200 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
     202              : 
     203         5876 :          DO ideriv = 1, nderivs
     204         4800 :             CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin), fm_struct)
     205              :          END DO
     206              :       END DO
     207              : 
     208              :       ! +++ allocate work matrices
     209         4428 :       ALLOCATE (S_mos_virt(nspins))
     210         2276 :       DO ispin = 1, nspins
     211         1200 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
     212         1200 :          CALL cp_fm_create(S_mos_virt(ispin), fm_struct)
     213              :          CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
     214              :                                       gs_mos(ispin)%mos_virt, &
     215              :                                       S_mos_virt(ispin), &
     216         2276 :                                       ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
     217              :       END DO
     218              : 
     219              :       ! check that the chosen dipole operator is consistent with the periodic boundary conditions used
     220         1076 :       CALL pw_env_get(pw_env, poisson_env=poisson_env)
     221         4304 :       ndim_periodic = COUNT(poisson_env%parameters%periodic == 1)
     222              : 
     223              :       ! select default for dipole form
     224         1076 :       IF (tddfpt_control%dipole_form == 0) THEN
     225          490 :          CALL get_qs_env(qs_env, dft_control=dft_control)
     226          490 :          IF (dft_control%qs_control%xtb) THEN
     227           32 :             IF (ndim_periodic == 0) THEN
     228            0 :                tddfpt_control%dipole_form = tddfpt_dipole_length
     229              :             ELSE
     230           32 :                tddfpt_control%dipole_form = tddfpt_dipole_berry
     231              :             END IF
     232              :          ELSE
     233          458 :             tddfpt_control%dipole_form = tddfpt_dipole_velocity
     234              :          END IF
     235              :       END IF
     236              : 
     237         1276 :       SELECT CASE (tddfpt_control%dipole_form)
     238              :       CASE (tddfpt_dipole_berry)
     239          200 :          IF (ndim_periodic /= 3) THEN
     240              :             CALL cp_warn(__LOCATION__, &
     241              :                          "Fully periodic Poisson solver (PERIODIC xyz) is needed "// &
     242            0 :                          "for oscillator strengths based on the Berry phase formula")
     243              :          END IF
     244              : 
     245          200 :          NULLIFY (berry_cossin_xyz)
     246              :          ! index: 1 = Re[exp(-i * G_t * t)],
     247              :          !        2 = Im[exp(-i * G_t * t)];
     248              :          ! t = x,y,z
     249          200 :          CALL dbcsr_allocate_matrix_set(berry_cossin_xyz, 2)
     250              : 
     251          600 :          DO i_cos_sin = 1, 2
     252          400 :             CALL dbcsr_init_p(berry_cossin_xyz(i_cos_sin)%matrix)
     253          600 :             CALL dbcsr_copy(berry_cossin_xyz(i_cos_sin)%matrix, matrix_s(1)%matrix)
     254              :          END DO
     255              : 
     256              :          ! +++ allocate berry-phase-related work matrices
     257         3000 :          ALLOCATE (gamma_00(nspins), gamma_inv_00(nspins), gamma_real_imag(2, nspins), opvec(2, nspins))
     258         1200 :          ALLOCATE (dBerry_mos_occ(nderivs, nspins))
     259          400 :          DO ispin = 1, nspins
     260          200 :             NULLIFY (fm_struct)
     261              :             CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_occ(ispin), &
     262          200 :                                      ncol_global=nmo_occ(ispin), context=blacs_env)
     263              : 
     264          200 :             CALL cp_cfm_create(gamma_00(ispin), fm_struct)
     265          200 :             CALL cp_cfm_create(gamma_inv_00(ispin), fm_struct)
     266              : 
     267          600 :             DO i_cos_sin = 1, 2
     268          600 :                CALL cp_fm_create(gamma_real_imag(i_cos_sin, ispin), fm_struct)
     269              :             END DO
     270          200 :             CALL cp_fm_struct_release(fm_struct)
     271              : 
     272              :             ! G_real C_0, G_imag C_0
     273          200 :             CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
     274          600 :             DO i_cos_sin = 1, 2
     275          600 :                CALL cp_fm_create(opvec(i_cos_sin, ispin), fm_struct)
     276              :             END DO
     277              : 
     278              :             ! dBerry * C_0
     279         1000 :             DO ideriv = 1, nderivs
     280          600 :                CALL cp_fm_create(dBerry_mos_occ(ideriv, ispin), fm_struct)
     281          800 :                CALL cp_fm_set_all(dBerry_mos_occ(ideriv, ispin), 0.0_dp)
     282              :             END DO
     283              :          END DO
     284              : 
     285          800 :          DO ideriv = 1, nderivs
     286         2400 :             kvec(:) = twopi*cell%h_inv(ideriv, :)
     287         1800 :             DO i_cos_sin = 1, 2
     288         1800 :                CALL dbcsr_set(berry_cossin_xyz(i_cos_sin)%matrix, 0.0_dp)
     289              :             END DO
     290              :             CALL build_berry_moment_matrix(qs_env, berry_cossin_xyz(1)%matrix, &
     291          600 :                                            berry_cossin_xyz(2)%matrix, kvec)
     292              : 
     293         1400 :             DO ispin = 1, nspins
     294              :                ! i_cos_sin = 1: cos (real) component; opvec(1) = gamma_real C_0
     295              :                ! i_cos_sin = 2: sin (imaginary) component; opvec(2) = gamma_imag C_0
     296         1800 :                DO i_cos_sin = 1, 2
     297              :                   CALL cp_dbcsr_sm_fm_multiply(berry_cossin_xyz(i_cos_sin)%matrix, &
     298              :                                                gs_mos(ispin)%mos_occ, &
     299              :                                                opvec(i_cos_sin, ispin), &
     300         1800 :                                                ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
     301              :                END DO
     302              : 
     303              :                CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, &
     304              :                                   1.0_dp, gs_mos(ispin)%mos_occ, opvec(1, ispin), &
     305          600 :                                   0.0_dp, gamma_real_imag(1, ispin))
     306              : 
     307              :                CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_occ(ispin), nao, &
     308              :                                   -1.0_dp, gs_mos(ispin)%mos_occ, opvec(2, ispin), &
     309          600 :                                   0.0_dp, gamma_real_imag(2, ispin))
     310              : 
     311              :                CALL cp_fm_to_cfm(msourcer=gamma_real_imag(1, ispin), &
     312              :                                  msourcei=gamma_real_imag(2, ispin), &
     313          600 :                                  mtarget=gamma_00(ispin))
     314              : 
     315              :                ! gamma_inv_00 = Q = [C_0^T (gamma_real - i gamma_imag) C_0] ^ {-1}
     316          600 :                CALL cp_cfm_set_all(gamma_inv_00(ispin), z_zero, z_one)
     317          600 :                CALL cp_cfm_solve(gamma_00(ispin), gamma_inv_00(ispin))
     318              : 
     319              :                CALL cp_cfm_to_fm(msource=gamma_inv_00(ispin), &
     320              :                                  mtargetr=gamma_real_imag(1, ispin), &
     321          600 :                                  mtargeti=gamma_real_imag(2, ispin))
     322              : 
     323              :                ! dBerry_mos_occ is identical to dBerry_psi0 from qs_linres_op % polar_operators()
     324              :                CALL parallel_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), &
     325              :                                   1.0_dp, opvec(1, ispin), gamma_real_imag(2, ispin), &
     326          600 :                                   0.0_dp, dipole_op_mos_occ(1, ispin))
     327              :                CALL parallel_gemm("N", "N", nao, nmo_occ(ispin), nmo_occ(ispin), &
     328              :                                   -1.0_dp, opvec(2, ispin), gamma_real_imag(1, ispin), &
     329          600 :                                   1.0_dp, dipole_op_mos_occ(1, ispin))
     330              : 
     331         3000 :                DO jderiv = 1, nderivs
     332              :                   CALL cp_fm_scale_and_add(1.0_dp, dBerry_mos_occ(jderiv, ispin), &
     333         2400 :                                            cell%hmat(jderiv, ideriv), dipole_op_mos_occ(1, ispin))
     334              :                END DO
     335              :             END DO
     336              :          END DO
     337              : 
     338              :          ! --- release berry-phase-related work matrices
     339          200 :          CALL cp_fm_release(opvec)
     340          200 :          CALL cp_fm_release(gamma_real_imag)
     341          400 :          DO ispin = nspins, 1, -1
     342          200 :             CALL cp_cfm_release(gamma_inv_00(ispin))
     343          400 :             CALL cp_cfm_release(gamma_00(ispin))
     344              :          END DO
     345          200 :          DEALLOCATE (gamma_00, gamma_inv_00)
     346          200 :          CALL dbcsr_deallocate_matrix_set(berry_cossin_xyz)
     347              : 
     348          200 :          NULLIFY (fm_struct)
     349          200 :          CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
     350          200 :          CALL cp_fm_create(wfm_ao_ao, fm_struct)
     351          200 :          CALL cp_fm_struct_release(fm_struct)
     352              : 
     353              :          ! trans_dipole = 2|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00) +
     354              :          !                2|e|/|G_mu| * Tr Imag(C_0^T * (gamma_real - i gamma_imag) * evects * gamma_inv_00) ,
     355              :          !
     356              :          ! Taking into account the symmetry of the matrices 'gamma_real' and 'gamma_imag' and the fact
     357              :          ! that the response wave-function is a real-valued function, the above expression can be simplified as
     358              :          ! trans_dipole = 4|e|/|G_mu| * Tr Imag(evects^T * (gamma_real - i gamma_imag) * C_0 * gamma_inv_00)
     359              :          !
     360              :          ! 1/|G_mu| = |lattice_vector_mu| / (2*pi) .
     361          400 :          DO ispin = 1, nspins
     362              :             ! wfm_ao_ao = S * mos_virt * mos_virt^T
     363              :             CALL parallel_gemm('N', 'T', nao, nao, nmo_virt(ispin), &
     364              :                                1.0_dp/twopi, S_mos_virt(ispin), gs_mos(ispin)%mos_virt, &
     365          200 :                                0.0_dp, wfm_ao_ao)
     366              : 
     367         1000 :             DO ideriv = 1, nderivs
     368              :                CALL parallel_gemm('N', 'N', nao, nmo_occ(ispin), nao, &
     369              :                                   1.0_dp, wfm_ao_ao, dBerry_mos_occ(ideriv, ispin), &
     370          800 :                                   0.0_dp, dipole_op_mos_occ(ideriv, ispin))
     371              :             END DO
     372              :          END DO
     373              : 
     374          200 :          CALL cp_fm_release(wfm_ao_ao)
     375          400 :          CALL cp_fm_release(dBerry_mos_occ)
     376              : 
     377              :       CASE (tddfpt_dipole_length)
     378            8 :          IF (ndim_periodic /= 0) THEN
     379              :             CALL cp_warn(__LOCATION__, &
     380              :                          "Non-periodic Poisson solver (PERIODIC none) is needed "// &
     381            0 :                          "for oscillator strengths based on the length operator")
     382              :          END IF
     383              : 
     384              :          ! compute components of the dipole operator in the length form
     385            8 :          NULLIFY (rRc_xyz)
     386            8 :          CALL dbcsr_allocate_matrix_set(rRc_xyz, nderivs)
     387              : 
     388           32 :          DO ideriv = 1, nderivs
     389           24 :             CALL dbcsr_init_p(rRc_xyz(ideriv)%matrix)
     390           32 :             CALL dbcsr_copy(rRc_xyz(ideriv)%matrix, matrix_s(1)%matrix)
     391              :          END DO
     392              : 
     393              :          CALL get_reference_point(reference_point, qs_env=qs_env, &
     394              :                                   reference=tddfpt_control%dipole_reference, &
     395            8 :                                   ref_point=tddfpt_control%dipole_ref_point)
     396              : 
     397              :          CALL rRc_xyz_ao(op=rRc_xyz, qs_env=qs_env, rc=reference_point, order=1, &
     398            8 :                          minimum_image=.FALSE., soft=.FALSE.)
     399              : 
     400            8 :          NULLIFY (fm_struct)
     401            8 :          CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
     402            8 :          CALL cp_fm_create(wfm_ao_ao, fm_struct)
     403            8 :          CALL cp_fm_struct_release(fm_struct)
     404              : 
     405           16 :          DO ispin = 1, nspins
     406            8 :             CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct)
     407            8 :             CALL cp_fm_create(rRc_mos_occ, fm_struct)
     408              : 
     409              :             ! wfm_ao_ao = S * mos_virt * mos_virt^T
     410              :             CALL parallel_gemm('N', 'T', nao, nao, nmo_virt(ispin), &
     411              :                                1.0_dp, S_mos_virt(ispin), gs_mos(ispin)%mos_virt, &
     412            8 :                                0.0_dp, wfm_ao_ao)
     413              : 
     414           32 :             DO ideriv = 1, nderivs
     415              :                CALL cp_dbcsr_sm_fm_multiply(rRc_xyz(ideriv)%matrix, &
     416              :                                             gs_mos(ispin)%mos_occ, &
     417              :                                             rRc_mos_occ, &
     418           24 :                                             ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
     419              : 
     420              :                CALL parallel_gemm('N', 'N', nao, nmo_occ(ispin), nao, &
     421              :                                   1.0_dp, wfm_ao_ao, rRc_mos_occ, &
     422           32 :                                   0.0_dp, dipole_op_mos_occ(ideriv, ispin))
     423              :             END DO
     424              : 
     425           24 :             CALL cp_fm_release(rRc_mos_occ)
     426              :          END DO
     427              : 
     428            8 :          CALL cp_fm_release(wfm_ao_ao)
     429            8 :          CALL dbcsr_deallocate_matrix_set(rRc_xyz)
     430              : 
     431              :       CASE (tddfpt_dipole_velocity)
     432              :          ! generate overlap derivatives
     433          868 :          CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb)
     434          868 :          NULLIFY (scrm)
     435              :          CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
     436              :                                    basis_type_a="ORB", basis_type_b="ORB", &
     437          868 :                                    sab_nl=sab_orb)
     438              : 
     439         1860 :          DO ispin = 1, nspins
     440          992 :             NULLIFY (fm_struct)
     441              :             CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), &
     442          992 :                                      ncol_global=nmo_occ(ispin), context=blacs_env)
     443          992 :             CALL cp_fm_create(ediff_inv, fm_struct)
     444          992 :             CALL cp_fm_create(wfm_mo_virt_mo_occ, fm_struct)
     445          992 :             CALL cp_fm_struct_release(fm_struct)
     446              : 
     447              :             CALL cp_fm_get_info(ediff_inv, nrow_local=nrows_local, ncol_local=ncols_local, &
     448          992 :                                 row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
     449          992 :             CALL cp_fm_get_info(wfm_mo_virt_mo_occ, local_data=local_data_wfm)
     450              : 
     451              : !$OMP       PARALLEL DO DEFAULT(NONE), &
     452              : !$OMP                PRIVATE(eval_occ, icol, irow), &
     453          992 : !$OMP                SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices)
     454              :             DO icol = 1, ncols_local
     455              :                ! E_occ_i ; imo_occ = col_indices(icol)
     456              :                eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol))
     457              : 
     458              :                DO irow = 1, nrows_local
     459              :                   ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i)
     460              :                   ! imo_virt = row_indices(irow)
     461              :                   local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ)
     462              :                END DO
     463              :             END DO
     464              : !$OMP       END PARALLEL DO
     465              : 
     466         3968 :             DO ideriv = 1, nderivs
     467              :                CALL cp_dbcsr_sm_fm_multiply(scrm(ideriv + 1)%matrix, &
     468              :                                             gs_mos(ispin)%mos_occ, &
     469              :                                             dipole_op_mos_occ(ideriv, ispin), &
     470         2976 :                                             ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
     471              : 
     472              :                CALL parallel_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, &
     473              :                                   1.0_dp, gs_mos(ispin)%mos_virt, dipole_op_mos_occ(ideriv, ispin), &
     474         2976 :                                   0.0_dp, wfm_mo_virt_mo_occ)
     475              : 
     476              :                ! in-place element-wise (Schur) product;
     477              :                ! avoid allocation of a temporary [nmo_virt x nmo_occ] matrix which is needed
     478              :                ! for cp_fm_schur_product() subroutine call
     479              : 
     480              : !$OMP          PARALLEL DO DEFAULT(NONE), &
     481              : !$OMP                   PRIVATE(icol, irow), &
     482         2976 : !$OMP                   SHARED(ispin, local_data_ediff, local_data_wfm, ncols_local, nrows_local)
     483              :                DO icol = 1, ncols_local
     484              :                   DO irow = 1, nrows_local
     485              :                      local_data_wfm(irow, icol) = local_data_wfm(irow, icol)*local_data_ediff(irow, icol)
     486              :                   END DO
     487              :                END DO
     488              : !$OMP          END PARALLEL DO
     489              : 
     490              :                CALL parallel_gemm('N', 'N', nao, nmo_occ(ispin), nmo_virt(ispin), &
     491              :                                   1.0_dp, S_mos_virt(ispin), wfm_mo_virt_mo_occ, &
     492         3968 :                                   0.0_dp, dipole_op_mos_occ(ideriv, ispin))
     493              :             END DO
     494              : 
     495          992 :             CALL cp_fm_release(wfm_mo_virt_mo_occ)
     496         3844 :             CALL cp_fm_release(ediff_inv)
     497              :          END DO
     498          868 :          CALL dbcsr_deallocate_matrix_set(scrm)
     499              : 
     500              :       CASE DEFAULT
     501         1084 :          CPABORT("Unimplemented form of the dipole operator")
     502              :       END SELECT
     503              : 
     504              :       ! --- release work matrices
     505         1076 :       CALL cp_fm_release(S_mos_virt)
     506              : 
     507         1076 :       CALL timestop(handle)
     508         3228 :    END SUBROUTINE tddfpt_dipole_operator
     509              : 
     510              : ! **************************************************************************************************
     511              : !> \brief Print final TDDFPT excitation energies and oscillator strengths.
     512              : !> \param log_unit           output unit
     513              : !> \param evects             TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
     514              : !>                           SIZE(evects,2) -- number of excited states to print)
     515              : !> \param evals              TDDFPT eigenvalues
     516              : !> \param ostrength          TDDFPT oscillator strength
     517              : !> \param mult               multiplicity
     518              : !> \param dipole_op_mos_occ  action of the dipole operator on the ground state wave function
     519              : !>                           [x,y,z ; spin]
     520              : !> \param dipole_form ...
     521              : !> \par History
     522              : !>    * 05.2016 created [Sergey Chulkov]
     523              : !>    * 06.2016 transition dipole moments and oscillator strengths [Sergey Chulkov]
     524              : !>    * 07.2016 spin-unpolarised electron density [Sergey Chulkov]
     525              : !>    * 08.2018 compute 'dipole_op_mos_occ' in a separate subroutine [Sergey Chulkov]
     526              : !> \note \parblock
     527              : !>       Adapted version of the subroutine find_contributions() which was originally created
     528              : !>       by Thomas Chassaing on 02.2005.
     529              : !>
     530              : !>       Transition dipole moment along direction 'd' is computed as following:
     531              : !>       \f[ t_d(spin) = Tr[evects^T dipole\_op\_mos\_occ(d, spin)] .\f]
     532              : !>       \endparblock
     533              : ! **************************************************************************************************
     534         2168 :    SUBROUTINE tddfpt_print_summary(log_unit, evects, evals, ostrength, mult, &
     535         1084 :                                    dipole_op_mos_occ, dipole_form)
     536              :       INTEGER, INTENT(in)                                :: log_unit
     537              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: evects
     538              :       REAL(kind=dp), DIMENSION(:), INTENT(in)            :: evals
     539              :       REAL(kind=dp), DIMENSION(:), INTENT(inout)         :: ostrength
     540              :       INTEGER, INTENT(in)                                :: mult
     541              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: dipole_op_mos_occ
     542              :       INTEGER, INTENT(in)                                :: dipole_form
     543              : 
     544              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_summary'
     545              : 
     546              :       CHARACTER(len=1)                                   :: lsd_str
     547              :       CHARACTER(len=20)                                  :: mult_str
     548              :       INTEGER                                            :: handle, ideriv, ispin, istate, nspins, &
     549              :                                                             nstates
     550              :       REAL(kind=dp)                                      :: osc_strength
     551         1084 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: trans_dipoles
     552              : 
     553         1084 :       CALL timeset(routineN, handle)
     554              : 
     555         1084 :       nspins = SIZE(evects, 1)
     556         1084 :       nstates = SIZE(evects, 2)
     557              : 
     558         1084 :       IF (nspins > 1) THEN
     559          124 :          lsd_str = 'U'
     560              :       ELSE
     561          960 :          lsd_str = 'R'
     562              :       END IF
     563              : 
     564              :       ! *** summary header ***
     565         1084 :       IF (log_unit > 0) THEN
     566          542 :          CALL integer_to_string(mult, mult_str)
     567          542 :          WRITE (log_unit, '(/,1X,A1,A,1X,A)') lsd_str, "-TDDFPT states of multiplicity", TRIM(mult_str)
     568          642 :          SELECT CASE (dipole_form)
     569              :          CASE (tddfpt_dipole_berry)
     570          100 :             WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using Berry operator formulation"
     571              :          CASE (tddfpt_dipole_length)
     572            7 :             WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using length formulation"
     573              :          CASE (tddfpt_dipole_velocity)
     574          435 :             WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using velocity formulation"
     575              :          CASE DEFAULT
     576          542 :             CPABORT("Unimplemented form of the dipole operator")
     577              :          END SELECT
     578              : 
     579          542 :          WRITE (log_unit, '(T10,A,T19,A,T37,A,T69,A)') "State", "Excitation", &
     580         1084 :             "Transition dipole (a.u.)", "Oscillator"
     581          542 :          WRITE (log_unit, '(T10,A,T19,A,T37,A,T49,A,T61,A,T67,A)') "number", "energy (eV)", &
     582         1084 :             "x", "y", "z", "strength (a.u.)"
     583          542 :          WRITE (log_unit, '(T10,72("-"))')
     584              :       END IF
     585              : 
     586              :       ! transition dipole moment
     587         4336 :       ALLOCATE (trans_dipoles(nstates, nderivs, nspins))
     588        15852 :       trans_dipoles(:, :, :) = 0.0_dp
     589              : 
     590              :       ! nspins == 1 .AND. mult == 3 : spin-flip transitions are forbidden due to symmetry reasons
     591         1084 :       IF (nspins > 1 .OR. mult == 1) THEN
     592         1940 :          DO ispin = 1, nspins
     593         5036 :             DO ideriv = 1, nderivs
     594              :                CALL cp_fm_trace(evects(ispin, :), dipole_op_mos_occ(ideriv, ispin), &
     595         4128 :                                 trans_dipoles(:, ideriv, ispin))
     596              :             END DO
     597              :          END DO
     598              : 
     599          908 :          IF (nspins == 1) THEN
     600         8974 :             trans_dipoles(:, :, 1) = SQRT(2.0_dp)*trans_dipoles(:, :, 1)
     601              :          ELSE
     602         1630 :             trans_dipoles(:, :, 1) = SQRT(trans_dipoles(:, :, 1)**2 + trans_dipoles(:, :, 2)**2)
     603              :          END IF
     604              :       END IF
     605              : 
     606              :       ! *** summary information ***
     607         4018 :       DO istate = 1, nstates
     608              :          osc_strength = 2.0_dp/3.0_dp*evals(istate)* &
     609         2934 :                         accurate_dot_product(trans_dipoles(istate, :, 1), trans_dipoles(istate, :, 1))
     610         2934 :          ostrength(istate) = osc_strength
     611         4018 :          IF (log_unit > 0) THEN
     612              :             WRITE (log_unit, '(1X,A,T9,I7,T19,F11.5,T31,3(1X,ES11.4E2),T69,ES12.5E2)') &
     613         1467 :                "TDDFPT|", istate, evals(istate)*evolt, trans_dipoles(istate, 1:nderivs, 1), osc_strength
     614              :          END IF
     615              :       END DO
     616              : 
     617              :       ! punch a checksum for the regs
     618         1084 :       IF (log_unit > 0) THEN
     619         2009 :          WRITE (log_unit, '(/,T2,A,E14.6)') 'TDDFPT : CheckSum  =', SQRT(SUM(evals**2))
     620              :       END IF
     621              : 
     622         1084 :       DEALLOCATE (trans_dipoles)
     623              : 
     624         1084 :       CALL timestop(handle)
     625         1084 :    END SUBROUTINE tddfpt_print_summary
     626              : 
     627              : ! **************************************************************************************************
     628              : !> \brief Print excitation analysis.
     629              : !> \param log_unit           output unit
     630              : !> \param evects             TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
     631              : !>                           SIZE(evects,2) -- number of excited states to print)
     632              : !> \param evals              TDDFPT eigenvalues
     633              : !> \param gs_mos             molecular orbitals optimised for the ground state
     634              : !> \param matrix_s           overlap matrix
     635              : !> \param min_amplitude      the smallest excitation amplitude to print
     636              : !> \par History
     637              : !>    * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
     638              : !>    * 08.2018 splited of from 'tddfpt_print_summary' [Sergey Chulkov]
     639              : ! **************************************************************************************************
     640         1084 :    SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, evals, gs_mos, matrix_s, min_amplitude)
     641              :       INTEGER, INTENT(in)                                :: log_unit
     642              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: evects
     643              :       REAL(kind=dp), DIMENSION(:), INTENT(in)            :: evals
     644              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     645              :          INTENT(in)                                      :: gs_mos
     646              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     647              :       REAL(kind=dp), INTENT(in)                          :: min_amplitude
     648              : 
     649              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_excitation_analysis'
     650              : 
     651              :       CHARACTER(len=5)                                   :: spin_label
     652              :       INTEGER                                            :: handle, icol, iproc, irow, ispin, &
     653              :                                                             istate, nao, ncols_local, nrows_local, &
     654              :                                                             nspins, nstates, state_spin
     655              :       INTEGER(kind=int_8)                                :: iexc, imo_occ, imo_virt, ind, nexcs, &
     656              :                                                             nexcs_local, nexcs_max_local, &
     657              :                                                             nmo_virt_occ, nmo_virt_occ_alpha
     658         1084 :       INTEGER(kind=int_8), ALLOCATABLE, DIMENSION(:)     :: inds_local, inds_recv, nexcs_recv
     659              :       INTEGER(kind=int_8), DIMENSION(1)                  :: nexcs_send
     660              :       INTEGER(kind=int_8), DIMENSION(maxspins)           :: nmo_occ8, nmo_virt8
     661         1084 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: inds
     662         1084 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     663              :       INTEGER, DIMENSION(maxspins)                       :: nmo_occ, nmo_virt
     664              :       LOGICAL                                            :: do_exc_analysis
     665         1084 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: weights_local, weights_neg_abs_recv, &
     666         1084 :                                                             weights_recv
     667              :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
     668         1084 :          POINTER                                         :: local_data
     669              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     670              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     671         1084 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: S_mos_virt, weights_fm
     672              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     673              :       TYPE(mp_request_type)                              :: send_handler, send_handler2
     674         1084 :       TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:)   :: recv_handlers, recv_handlers2
     675              : 
     676         1084 :       CALL timeset(routineN, handle)
     677              : 
     678         1084 :       nspins = SIZE(evects, 1)
     679         1084 :       nstates = SIZE(evects, 2)
     680         1084 :       do_exc_analysis = min_amplitude < 1.0_dp
     681              : 
     682         1084 :       CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env, para_env=para_env)
     683         1084 :       CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
     684              : 
     685         2292 :       DO ispin = 1, nspins
     686         1208 :          nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
     687         1208 :          nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
     688         1208 :          nmo_occ8(ispin) = SIZE(gs_mos(ispin)%evals_occ, kind=int_8)
     689         2292 :          nmo_virt8(ispin) = SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
     690              :       END DO
     691              : 
     692              :       ! *** excitation analysis ***
     693         1084 :       IF (do_exc_analysis) THEN
     694         1084 :          CPASSERT(log_unit <= 0 .OR. para_env%is_source())
     695         1084 :          nmo_virt_occ_alpha = INT(nmo_virt(1), int_8)*INT(nmo_occ(1), int_8)
     696              : 
     697         1084 :          IF (log_unit > 0) THEN
     698          542 :             WRITE (log_unit, "(1X,A)") "", &
     699          542 :                "-------------------------------------------------------------------------------", &
     700          542 :                "-                            Excitation analysis                              -", &
     701         1084 :                "-------------------------------------------------------------------------------"
     702          542 :             WRITE (log_unit, '(8X,A,T27,A,T49,A,T69,A)') "State", "Occupied", "Virtual", "Excitation"
     703          542 :             WRITE (log_unit, '(8X,A,T28,A,T49,A,T69,A)') "number", "orbital", "orbital", "amplitude"
     704          542 :             WRITE (log_unit, '(1X,79("-"))')
     705              : 
     706          542 :             IF (nspins == 1) THEN
     707          480 :                state_spin = 1
     708          480 :                spin_label = '     '
     709              :             END IF
     710              :          END IF
     711              : 
     712         6752 :          ALLOCATE (S_mos_virt(nspins), weights_fm(nspins))
     713         2292 :          DO ispin = 1, nspins
     714         1208 :             CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
     715         1208 :             CALL cp_fm_create(S_mos_virt(ispin), fm_struct)
     716              :             CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
     717              :                                          gs_mos(ispin)%mos_virt, &
     718              :                                          S_mos_virt(ispin), &
     719         1208 :                                          ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
     720              : 
     721         1208 :             NULLIFY (fm_struct)
     722              :             CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), &
     723         1208 :                                      context=blacs_env)
     724         1208 :             CALL cp_fm_create(weights_fm(ispin), fm_struct)
     725         2292 :             CALL cp_fm_struct_release(fm_struct)
     726              :          END DO
     727              : 
     728              :          nexcs_max_local = 0
     729         2292 :          DO ispin = 1, nspins
     730         1208 :             CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local)
     731         2292 :             nexcs_max_local = nexcs_max_local + INT(nrows_local, int_8)*INT(ncols_local, int_8)
     732              :          END DO
     733              : 
     734         4336 :          ALLOCATE (weights_local(nexcs_max_local), inds_local(nexcs_max_local))
     735              : 
     736         4018 :          DO istate = 1, nstates
     737              :             nexcs_local = 0
     738              :             nmo_virt_occ = 0
     739              : 
     740              :             ! analyse matrix elements locally and transfer only significant
     741              :             ! excitations to the master node for subsequent ordering
     742         6246 :             DO ispin = 1, nspins
     743              :                ! compute excitation amplitudes
     744              :                CALL parallel_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, 1.0_dp, S_mos_virt(ispin), &
     745         3312 :                                   evects(ispin, istate), 0.0_dp, weights_fm(ispin))
     746              : 
     747              :                CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
     748         3312 :                                    row_indices=row_indices, col_indices=col_indices, local_data=local_data)
     749              : 
     750              :                ! locate single excitations with significant amplitudes (>= min_amplitude)
     751        18376 :                DO icol = 1, ncols_local
     752       198322 :                   DO irow = 1, nrows_local
     753       195010 :                      IF (ABS(local_data(irow, icol)) >= min_amplitude) THEN
     754              :                         ! number of non-negligible excitations
     755         2303 :                         nexcs_local = nexcs_local + 1
     756              :                         ! excitation amplitude
     757         2303 :                         weights_local(nexcs_local) = local_data(irow, icol)
     758              :                         ! index of single excitation (ivirt, iocc, ispin) in compressed form
     759              :                         inds_local(nexcs_local) = nmo_virt_occ + INT(row_indices(irow), int_8) + &
     760         2303 :                                                   INT(col_indices(icol) - 1, int_8)*nmo_virt8(ispin)
     761              :                      END IF
     762              :                   END DO
     763              :                END DO
     764              : 
     765         9558 :                nmo_virt_occ = nmo_virt_occ + nmo_virt8(ispin)*nmo_occ8(ispin)
     766              :             END DO
     767              : 
     768         2934 :             IF (para_env%is_source()) THEN
     769              :                ! master node
     770        14670 :                ALLOCATE (nexcs_recv(para_env%num_pe), recv_handlers(para_env%num_pe), recv_handlers2(para_env%num_pe))
     771              : 
     772              :                ! collect number of non-negligible excitations from other nodes
     773         4401 :                DO iproc = 1, para_env%num_pe
     774         4401 :                   IF (iproc - 1 /= para_env%mepos) THEN
     775         1467 :                      CALL para_env%irecv(nexcs_recv(iproc:iproc), iproc - 1, recv_handlers(iproc), 0)
     776              :                   ELSE
     777         1467 :                      nexcs_recv(iproc) = nexcs_local
     778              :                   END IF
     779              :                END DO
     780              : 
     781         4401 :                DO iproc = 1, para_env%num_pe
     782         2934 :                   IF (iproc - 1 /= para_env%mepos) &
     783         2934 :                      CALL recv_handlers(iproc)%wait()
     784              :                END DO
     785              : 
     786              :                ! compute total number of non-negligible excitations
     787         1467 :                nexcs = 0
     788         4401 :                DO iproc = 1, para_env%num_pe
     789         4401 :                   nexcs = nexcs + nexcs_recv(iproc)
     790              :                END DO
     791              : 
     792              :                ! receive indices and amplitudes of selected excitations
     793         5868 :                ALLOCATE (weights_recv(nexcs), weights_neg_abs_recv(nexcs))
     794         5868 :                ALLOCATE (inds_recv(nexcs), inds(nexcs))
     795              : 
     796         4401 :                nmo_virt_occ = 0
     797         4401 :                DO iproc = 1, para_env%num_pe
     798         4401 :                   IF (nexcs_recv(iproc) > 0) THEN
     799         1538 :                      IF (iproc - 1 /= para_env%mepos) THEN
     800              :                         ! excitation amplitudes
     801              :                         CALL para_env%irecv(weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
     802          227 :                                             iproc - 1, recv_handlers(iproc), 1)
     803              :                         ! compressed indices
     804              :                         CALL para_env%irecv(inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)), &
     805          227 :                                             iproc - 1, recv_handlers2(iproc), 2)
     806              :                      ELSE
     807              :                         ! data on master node
     808         3350 :                         weights_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = weights_local(1:nexcs_recv(iproc))
     809         3350 :                         inds_recv(nmo_virt_occ + 1:nmo_virt_occ + nexcs_recv(iproc)) = inds_local(1:nexcs_recv(iproc))
     810              :                      END IF
     811              : 
     812         1538 :                      nmo_virt_occ = nmo_virt_occ + nexcs_recv(iproc)
     813              :                   END IF
     814              :                END DO
     815              : 
     816         4401 :                DO iproc = 1, para_env%num_pe
     817         4401 :                   IF (iproc - 1 /= para_env%mepos .AND. nexcs_recv(iproc) > 0) THEN
     818          227 :                      CALL recv_handlers(iproc)%wait()
     819          227 :                      CALL recv_handlers2(iproc)%wait()
     820              :                   END IF
     821              :                END DO
     822              : 
     823         1467 :                DEALLOCATE (nexcs_recv, recv_handlers, recv_handlers2)
     824              :             ELSE
     825              :                ! working node: send the number of selected excited states to the master node
     826         1467 :                nexcs_send(1) = nexcs_local
     827         1467 :                CALL para_env%isend(nexcs_send, para_env%source, send_handler, 0)
     828         1467 :                CALL send_handler%wait()
     829              : 
     830         1467 :                IF (nexcs_local > 0) THEN
     831              :                   ! send excitation amplitudes
     832          227 :                   CALL para_env%isend(weights_local(1:nexcs_local), para_env%source, send_handler, 1)
     833              :                   ! send compressed indices
     834          227 :                   CALL para_env%isend(inds_local(1:nexcs_local), para_env%source, send_handler2, 2)
     835              : 
     836          227 :                   CALL send_handler%wait()
     837          227 :                   CALL send_handler2%wait()
     838              :                END IF
     839              :             END IF
     840              : 
     841              :             ! sort non-negligible excitations on the master node according to their amplitudes,
     842              :             ! uncompress indices and print summary information
     843         2934 :             IF (para_env%is_source() .AND. log_unit > 0) THEN
     844         3770 :                weights_neg_abs_recv(:) = -ABS(weights_recv)
     845         1467 :                CALL sort(weights_neg_abs_recv, INT(nexcs), inds)
     846              : 
     847         1467 :                WRITE (log_unit, '(T7,I8,F10.5,A)') istate, evals(istate)*evolt, " eV"
     848              : 
     849         3770 :                DO iexc = 1, nexcs
     850         2303 :                   ind = inds_recv(inds(iexc)) - 1
     851         2303 :                   IF (nspins > 1) THEN
     852          390 :                      IF (ind < nmo_virt_occ_alpha) THEN
     853          176 :                         state_spin = 1
     854          176 :                         spin_label = '(alp)'
     855              :                      ELSE
     856          214 :                         state_spin = 2
     857          214 :                         ind = ind - nmo_virt_occ_alpha
     858          214 :                         spin_label = '(bet)'
     859              :                      END IF
     860              :                   END IF
     861              : 
     862         2303 :                   imo_occ = ind/nmo_virt8(state_spin) + 1
     863         2303 :                   imo_virt = MOD(ind, nmo_virt8(state_spin)) + 1
     864              : 
     865         2303 :                   WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
     866         6073 :                      nmo_occ8(state_spin) + imo_virt, spin_label, weights_recv(inds(iexc))
     867              :                END DO
     868              :             END IF
     869              : 
     870              :             ! deallocate temporary arrays
     871         2934 :             IF (para_env%is_source()) &
     872         2551 :                DEALLOCATE (weights_recv, weights_neg_abs_recv, inds_recv, inds)
     873              :          END DO
     874              : 
     875         1084 :          DEALLOCATE (weights_local, inds_local)
     876         1084 :          IF (log_unit > 0) THEN
     877              :             WRITE (log_unit, "(1X,A)") &
     878          542 :                "-------------------------------------------------------------------------------"
     879              :          END IF
     880              :       END IF
     881              : 
     882         1084 :       CALL cp_fm_release(weights_fm)
     883         1084 :       CALL cp_fm_release(S_mos_virt)
     884              : 
     885         1084 :       CALL timestop(handle)
     886              : 
     887         2168 :    END SUBROUTINE tddfpt_print_excitation_analysis
     888              : 
     889              : ! **************************************************************************************************
     890              : !> \brief Print natural transition orbital analysis.
     891              : !> \param qs_env             Information on Kinds and Particles
     892              : !> \param evects             TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
     893              : !>                           SIZE(evects,2) -- number of excited states to print)
     894              : !> \param evals              TDDFPT eigenvalues
     895              : !> \param ostrength ...
     896              : !> \param gs_mos             molecular orbitals optimised for the ground state
     897              : !> \param matrix_s           overlap matrix
     898              : !> \param print_section      ...
     899              : !> \par History
     900              : !>    * 06.2019 created [JGH]
     901              : ! **************************************************************************************************
     902         1084 :    SUBROUTINE tddfpt_print_nto_analysis(qs_env, evects, evals, ostrength, gs_mos, matrix_s, print_section)
     903              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     904              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: evects
     905              :       REAL(kind=dp), DIMENSION(:), INTENT(in)            :: evals, ostrength
     906              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
     907              :          INTENT(in)                                      :: gs_mos
     908              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     909              :       TYPE(section_vals_type), POINTER                   :: print_section
     910              : 
     911              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_nto_analysis'
     912              :       INTEGER, PARAMETER                                 :: ntomax = 10
     913              : 
     914              :       CHARACTER(LEN=20), DIMENSION(2)                    :: nto_name
     915              :       INTEGER                                            :: handle, i, ia, icg, iounit, ispin, &
     916              :                                                             istate, j, nao, nlist, nmax, nmo, &
     917              :                                                             nnto, nspins, nstates
     918              :       INTEGER, DIMENSION(2)                              :: iv
     919              :       INTEGER, DIMENSION(2, ntomax)                      :: ia_index
     920         1084 :       INTEGER, DIMENSION(:), POINTER                     :: slist, stride
     921              :       LOGICAL                                            :: append_cube, cube_file, explicit
     922              :       REAL(KIND=dp)                                      :: os_threshold, sume, threshold
     923         1084 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigvals
     924         1084 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: eigenvalues
     925              :       REAL(KIND=dp), DIMENSION(ntomax)                   :: ia_eval
     926              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_mo_struct, fm_struct
     927              :       TYPE(cp_fm_type)                                   :: Sev, smat, tmat, wmat, work, wvec
     928         1084 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: teig
     929              :       TYPE(cp_logger_type), POINTER                      :: logger
     930         1084 :       TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:)       :: nto_set
     931         1084 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     932         1084 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     933              :       TYPE(section_vals_type), POINTER                   :: molden_section, nto_section
     934              : 
     935         1084 :       CALL timeset(routineN, handle)
     936              : 
     937         1084 :       logger => cp_get_default_logger()
     938         1084 :       iounit = cp_logger_get_default_io_unit(logger)
     939              : 
     940         1084 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, &
     941              :                                            "NTO_ANALYSIS"), cp_p_file)) THEN
     942              : 
     943          144 :          CALL cite_reference(Martin2003)
     944              : 
     945          144 :          CALL section_vals_val_get(print_section, "NTO_ANALYSIS%THRESHOLD", r_val=threshold)
     946          144 :          CALL section_vals_val_get(print_section, "NTO_ANALYSIS%INTENSITY_THRESHOLD", r_val=os_threshold)
     947          144 :          CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", EXPLICIT=explicit)
     948              : 
     949          144 :          IF (explicit) THEN
     950            4 :             CALL section_vals_val_get(print_section, "NTO_ANALYSIS%STATE_LIST", i_vals=slist)
     951            4 :             nlist = SIZE(slist)
     952              :          ELSE
     953              :             nlist = 0
     954              :          END IF
     955              : 
     956          144 :          IF (iounit > 0) THEN
     957           72 :             WRITE (iounit, "(1X,A)") "", &
     958           72 :                "-------------------------------------------------------------------------------", &
     959           72 :                "-                            Natural Orbital analysis                         -", &
     960          144 :                "-------------------------------------------------------------------------------"
     961              :          END IF
     962              : 
     963          144 :          nspins = SIZE(evects, 1)
     964          144 :          nstates = SIZE(evects, 2)
     965          144 :          CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
     966              : 
     967          344 :          DO istate = 1, nstates
     968          200 :             IF (os_threshold > ostrength(istate)) THEN
     969           44 :                IF (iounit > 0) THEN
     970           22 :                   WRITE (iounit, "(1X,A,I6)") "  Skipping state ", istate
     971              :                END IF
     972              :                CYCLE
     973              :             END IF
     974          156 :             IF (nlist > 0) THEN
     975           12 :                IF (.NOT. ANY(slist == istate)) THEN
     976            0 :                   IF (iounit > 0) THEN
     977            0 :                      WRITE (iounit, "(1X,A,I6)") "  Skipping state ", istate
     978              :                   END IF
     979              :                   CYCLE
     980              :                END IF
     981              :             END IF
     982          156 :             IF (iounit > 0) THEN
     983           78 :                WRITE (iounit, "(1X,A,I6,T30,F10.5,A)") "  STATE NR. ", istate, evals(istate)*evolt, " eV"
     984              :             END IF
     985              :             nmax = 0
     986          326 :             DO ispin = 1, nspins
     987          170 :                CALL cp_fm_get_info(evects(ispin, istate), matrix_struct=fm_struct, ncol_global=nmo)
     988          326 :                nmax = MAX(nmax, nmo)
     989              :             END DO
     990          624 :             ALLOCATE (eigenvalues(nmax, nspins))
     991         1434 :             eigenvalues = 0.0_dp
     992              :             ! SET 1: Hole states
     993              :             ! SET 2: Particle states
     994          156 :             nto_name(1) = 'Hole_states'
     995          156 :             nto_name(2) = 'Particle_states'
     996          468 :             ALLOCATE (nto_set(2))
     997          468 :             DO i = 1, 2
     998          312 :                CALL allocate_mo_set(nto_set(i), nao, ntomax, 0, 0.0_dp, 1.0_dp, 0.0_dp)
     999          312 :                CALL cp_fm_get_info(evects(1, istate), matrix_struct=fm_struct)
    1000              :                CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
    1001          312 :                                         ncol_global=ntomax)
    1002          312 :                CALL cp_fm_create(tmat, fm_mo_struct)
    1003          312 :                CALL init_mo_set(nto_set(i), fm_ref=tmat, name=nto_name(i))
    1004          312 :                CALL cp_fm_release(tmat)
    1005          780 :                CALL cp_fm_struct_release(fm_mo_struct)
    1006              :             END DO
    1007              :             !
    1008          638 :             ALLOCATE (teig(nspins))
    1009              :             ! hole states
    1010              :             ! Diagonalize X(T)*S*X
    1011          326 :             DO ispin = 1, nspins
    1012              :                ASSOCIATE (ev => evects(ispin, istate))
    1013          170 :                   CALL cp_fm_get_info(ev, matrix_struct=fm_struct, ncol_global=nmo)
    1014          170 :                   CALL cp_fm_create(Sev, fm_struct)
    1015              :                   CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
    1016          170 :                                            nrow_global=nmo, ncol_global=nmo)
    1017          170 :                   CALL cp_fm_create(tmat, fm_mo_struct)
    1018          170 :                   CALL cp_fm_create(teig(ispin), fm_mo_struct)
    1019          170 :                   CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, Sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
    1020          170 :                   CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, ev, Sev, 0.0_dp, tmat)
    1021              :                END ASSOCIATE
    1022              : 
    1023          170 :                CALL choose_eigv_solver(tmat, teig(ispin), eigenvalues(1:nmo, ispin))
    1024              : 
    1025          170 :                CALL cp_fm_struct_release(fm_mo_struct)
    1026          170 :                CALL cp_fm_release(tmat)
    1027          666 :                CALL cp_fm_release(Sev)
    1028              :             END DO
    1029              :             ! find major determinants i->a
    1030          156 :             ia_index = 0
    1031          156 :             sume = 0.0_dp
    1032          156 :             nnto = 0
    1033          186 :             DO i = 1, ntomax
    1034          186 :                iv = MAXLOC(eigenvalues)
    1035          186 :                ia_eval(i) = eigenvalues(iv(1), iv(2))
    1036          558 :                ia_index(1:2, i) = iv(1:2)
    1037          186 :                sume = sume + ia_eval(i)
    1038          186 :                eigenvalues(iv(1), iv(2)) = 0.0_dp
    1039          186 :                nnto = nnto + 1
    1040          186 :                IF (sume > threshold) EXIT
    1041              :             END DO
    1042              :             ! store hole states
    1043          156 :             CALL set_mo_set(nto_set(1), nmo=nnto)
    1044          342 :             DO i = 1, nnto
    1045          186 :                ia = ia_index(1, i)
    1046          186 :                ispin = ia_index(2, i)
    1047          186 :                CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=nmo)
    1048          186 :                CALL cp_fm_get_info(teig(ispin), matrix_struct=fm_struct)
    1049              :                CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
    1050          186 :                                         nrow_global=nmo, ncol_global=1)
    1051          186 :                CALL cp_fm_create(tmat, fm_mo_struct)
    1052          186 :                CALL cp_fm_struct_release(fm_mo_struct)
    1053          186 :                CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=fm_struct)
    1054              :                CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
    1055          186 :                                         ncol_global=1)
    1056          186 :                CALL cp_fm_create(wvec, fm_mo_struct)
    1057          186 :                CALL cp_fm_struct_release(fm_mo_struct)
    1058          186 :                CALL cp_fm_to_fm(teig(ispin), tmat, 1, ia, 1)
    1059              :                CALL parallel_gemm('N', 'N', nao, 1, nmo, 1.0_dp, gs_mos(ispin)%mos_occ, &
    1060          186 :                                   tmat, 0.0_dp, wvec)
    1061          186 :                CALL cp_fm_to_fm(wvec, nto_set(1)%mo_coeff, 1, 1, i)
    1062          186 :                CALL cp_fm_release(wvec)
    1063          900 :                CALL cp_fm_release(tmat)
    1064              :             END DO
    1065              :             ! particle states
    1066              :             ! Solve generalized eigenvalue equation:  (S*X)*(S*X)(T)*v = lambda*S*v
    1067          156 :             CALL set_mo_set(nto_set(2), nmo=nnto)
    1068          326 :             DO ispin = 1, nspins
    1069              :                ASSOCIATE (ev => evects(ispin, istate))
    1070          170 :                   CALL cp_fm_get_info(ev, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
    1071          510 :                   ALLOCATE (eigvals(nao))
    1072         4534 :                   eigvals = 0.0_dp
    1073          170 :                   CALL cp_fm_create(Sev, fm_struct)
    1074          340 :                   CALL cp_dbcsr_sm_fm_multiply(matrix_s, ev, Sev, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
    1075              :                END ASSOCIATE
    1076              :                CALL cp_fm_struct_create(fmstruct=fm_mo_struct, template_fmstruct=fm_struct, &
    1077          170 :                                         nrow_global=nao, ncol_global=nao)
    1078          170 :                CALL cp_fm_create(tmat, fm_mo_struct)
    1079          170 :                CALL cp_fm_create(smat, fm_mo_struct)
    1080          170 :                CALL cp_fm_create(wmat, fm_mo_struct)
    1081          170 :                CALL cp_fm_create(work, fm_mo_struct)
    1082          170 :                CALL cp_fm_struct_release(fm_mo_struct)
    1083          170 :                CALL copy_dbcsr_to_fm(matrix_s, smat)
    1084          170 :                CALL parallel_gemm('N', 'T', nao, nao, nmo, 1.0_dp, Sev, Sev, 0.0_dp, tmat)
    1085          170 :                CALL cp_fm_geeig(tmat, smat, wmat, eigvals, work)
    1086          376 :                DO i = 1, nnto
    1087          376 :                   IF (ispin == ia_index(2, i)) THEN
    1088          186 :                      icg = 0
    1089         4584 :                      DO j = 1, nao
    1090         4584 :                         IF (ABS(eigvals(j) - ia_eval(i)) < 1.E-6_dp) THEN
    1091          186 :                            icg = j
    1092          186 :                            EXIT
    1093              :                         END IF
    1094              :                      END DO
    1095          186 :                      IF (icg == 0) THEN
    1096              :                         CALL cp_warn(__LOCATION__, &
    1097            0 :                                      "Could not locate particle state associated with hole state.")
    1098              :                      ELSE
    1099          186 :                         CALL cp_fm_to_fm(wmat, nto_set(2)%mo_coeff, 1, icg, i)
    1100              :                      END IF
    1101              :                   END IF
    1102              :                END DO
    1103          170 :                DEALLOCATE (eigvals)
    1104          170 :                CALL cp_fm_release(Sev)
    1105          170 :                CALL cp_fm_release(tmat)
    1106          170 :                CALL cp_fm_release(smat)
    1107          170 :                CALL cp_fm_release(wmat)
    1108          496 :                CALL cp_fm_release(work)
    1109              :             END DO
    1110              :             ! print
    1111          156 :             IF (iounit > 0) THEN
    1112           78 :                sume = 0.0_dp
    1113          171 :                DO i = 1, nnto
    1114           93 :                   sume = sume + ia_eval(i)
    1115              :                   WRITE (iounit, "(T6,A,i2,T30,A,i1,T42,A,F8.5,T63,A,F8.5)") &
    1116           93 :                      "Particle-Hole state:", i, " Spin:", ia_index(2, i), &
    1117          264 :                      "Eigenvalue:", ia_eval(i), " Sum Eigv:", sume
    1118              :                END DO
    1119              :             END IF
    1120              :             ! Cube and Molden files
    1121          156 :             nto_section => section_vals_get_subs_vals(print_section, "NTO_ANALYSIS")
    1122          156 :             CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
    1123          156 :             CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
    1124          156 :             CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
    1125          156 :             IF (cube_file) THEN
    1126           16 :                CALL print_nto_cubes(qs_env, nto_set, istate, stride, append_cube, nto_section)
    1127              :             END IF
    1128          156 :             CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set)
    1129          156 :             molden_section => section_vals_get_subs_vals(print_section, "MOS_MOLDEN")
    1130          156 :             CALL write_mos_molden(nto_set, qs_kind_set, particle_set, molden_section)
    1131              :             !
    1132          156 :             DEALLOCATE (eigenvalues)
    1133          156 :             CALL cp_fm_release(teig)
    1134              :             !
    1135          468 :             DO i = 1, 2
    1136          468 :                CALL deallocate_mo_set(nto_set(i))
    1137              :             END DO
    1138          612 :             DEALLOCATE (nto_set)
    1139              :          END DO
    1140              : 
    1141          144 :          IF (iounit > 0) THEN
    1142              :             WRITE (iounit, "(1X,A)") &
    1143           72 :                "-------------------------------------------------------------------------------"
    1144              :          END IF
    1145              : 
    1146              :       END IF
    1147              : 
    1148         1084 :       CALL timestop(handle)
    1149              : 
    1150         2168 :    END SUBROUTINE tddfpt_print_nto_analysis
    1151              : 
    1152              : ! **************************************************************************************************
    1153              : !> \brief Print exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
    1154              : !> \param log_unit                              output unit
    1155              : !> \param evects                                TDDFPT trial vectors (SIZE(evects,1) -- number of spins;
    1156              : !>                                              SIZE(evects,2) -- number of excited states to print)
    1157              : !> \param gs_mos                                molecular orbitals optimised for the ground state
    1158              : !> \param matrix_s                              overlap matrix
    1159              : !> \param do_directional_exciton_descriptors    flag for computing descriptors for each (cartesian) direction
    1160              : !> \param qs_env                                Information on particles/geometry
    1161              : !> \par History
    1162              : !>    * 12.2024 created as 'tddfpt_print_exciton_descriptors' [Maximilian Graml]
    1163              : ! **************************************************************************************************
    1164            2 :    SUBROUTINE tddfpt_print_exciton_descriptors(log_unit, evects, gs_mos, matrix_s, &
    1165              :                                                do_directional_exciton_descriptors, qs_env)
    1166              :       INTEGER, INTENT(in)                                :: log_unit
    1167              :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in)      :: evects
    1168              :       TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
    1169              :          INTENT(in)                                      :: gs_mos
    1170              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
    1171              :       LOGICAL, INTENT(IN) :: do_directional_exciton_descriptors
    1172              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
    1173              : 
    1174              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_exciton_descriptors'
    1175              : 
    1176              :       CHARACTER(LEN=4)                                   :: prefix_output
    1177              :       INTEGER                                            :: handle, ispin, istate, n_moments_quad, &
    1178              :                                                             nao, nspins, nstates
    1179              :       INTEGER, DIMENSION(maxspins)                       :: nmo_occ, nmo_virt
    1180              :       LOGICAL                                            :: print_checkvalue
    1181            2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ref_point_multipole
    1182              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1183              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_mo_coeff, &
    1184              :                                                             fm_struct_S_mos_virt, fm_struct_X_ia_n
    1185            2 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: eigvec_X_ia_n, fm_multipole_ab, &
    1186            2 :                                                             fm_multipole_ai, fm_multipole_ij, &
    1187            2 :                                                             S_mos_virt
    1188            2 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mo_coeff
    1189              :       TYPE(exciton_descr_type), ALLOCATABLE, &
    1190            2 :          DIMENSION(:)                                    :: exc_descr
    1191              : 
    1192            2 :       CALL timeset(routineN, handle)
    1193              : 
    1194            2 :       nspins = SIZE(evects, 1)
    1195            2 :       nstates = SIZE(evects, 2)
    1196              : 
    1197            2 :       CPASSERT(nspins == 1) ! Other spins are not yet implemented for exciton descriptors
    1198              : 
    1199            2 :       CALL cp_fm_get_info(gs_mos(1)%mos_occ, context=blacs_env)
    1200            2 :       CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
    1201              : 
    1202            4 :       DO ispin = 1, nspins
    1203            2 :          nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
    1204            4 :          nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
    1205              :       END DO
    1206              : 
    1207              :       ! Prepare fm with all MO coefficents, i.e. nao x nao
    1208            8 :       ALLOCATE (mo_coeff(nspins))
    1209              :       CALL cp_fm_struct_create(fm_struct_mo_coeff, nrow_global=nao, ncol_global=nao, &
    1210            2 :                                context=blacs_env)
    1211            4 :       DO ispin = 1, nspins
    1212            2 :          CALL cp_fm_create(mo_coeff(ispin), fm_struct_mo_coeff)
    1213              :          CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_occ, &
    1214              :                                          mo_coeff(ispin), &
    1215              :                                          nao, &
    1216              :                                          nmo_occ(ispin), &
    1217              :                                          1, &
    1218              :                                          1, &
    1219              :                                          1, &
    1220              :                                          1, &
    1221            2 :                                          blacs_env)
    1222              :          CALL cp_fm_to_fm_submat_general(gs_mos(ispin)%mos_virt, &
    1223              :                                          mo_coeff(ispin), &
    1224              :                                          nao, &
    1225              :                                          nmo_virt(ispin), &
    1226              :                                          1, &
    1227              :                                          1, &
    1228              :                                          1, &
    1229              :                                          nmo_occ(ispin) + 1, &
    1230            4 :                                          blacs_env)
    1231              :       END DO
    1232            2 :       CALL cp_fm_struct_release(fm_struct_mo_coeff)
    1233              : 
    1234              :       ! Compute multipole moments
    1235              :       ! fm_multipole_XY have structure inherited by libint, i.e. x, y, z, xx, xy, xz, yy, yz, zz
    1236            2 :       n_moments_quad = 9
    1237            2 :       ALLOCATE (ref_point_multipole(3))
    1238           20 :       ALLOCATE (fm_multipole_ij(n_moments_quad))
    1239           20 :       ALLOCATE (fm_multipole_ab(n_moments_quad))
    1240           20 :       ALLOCATE (fm_multipole_ai(n_moments_quad))
    1241              : 
    1242              :       CALL get_multipoles_mo(fm_multipole_ai, fm_multipole_ij, fm_multipole_ab, &
    1243              :                              qs_env, mo_coeff, ref_point_multipole, 2, &
    1244            2 :                              nmo_occ(1), nmo_virt(1), blacs_env)
    1245              : 
    1246            2 :       CALL cp_fm_release(mo_coeff)
    1247              : 
    1248              :       ! Compute eigenvector X of the Casida equation from trial vectors
    1249           14 :       ALLOCATE (S_mos_virt(nspins), eigvec_X_ia_n(nspins))
    1250            4 :       DO ispin = 1, nspins
    1251            2 :          CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct_S_mos_virt)
    1252            2 :          CALL cp_fm_create(S_mos_virt(ispin), fm_struct_S_mos_virt)
    1253            2 :          NULLIFY (fm_struct_S_mos_virt)
    1254              :          CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
    1255              :                                       gs_mos(ispin)%mos_virt, &
    1256              :                                       S_mos_virt(ispin), &
    1257            2 :                                       ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
    1258              : 
    1259              :          CALL cp_fm_struct_create(fm_struct_X_ia_n, nrow_global=nmo_occ(ispin), ncol_global=nmo_virt(ispin), &
    1260            2 :                                   context=blacs_env)
    1261            2 :          CALL cp_fm_create(eigvec_X_ia_n(ispin), fm_struct_X_ia_n)
    1262            4 :          CALL cp_fm_struct_release(fm_struct_X_ia_n)
    1263              :       END DO
    1264           78 :       ALLOCATE (exc_descr(nstates))
    1265           12 :       DO istate = 1, nstates
    1266           22 :          DO ispin = 1, nspins
    1267           10 :             CALL cp_fm_set_all(eigvec_X_ia_n(ispin), 0.0_dp)
    1268              :             ! compute eigenvectors X of the TDA equation
    1269              :             ! Reshuffle multiplication from
    1270              :             ! X_ai = S_µa ^T * C_µi
    1271              :             ! to
    1272              :             ! X_ia = C_µi ^T * S_µa
    1273              :             ! for compatibility with the structure needed for get_exciton_descriptors of bse_properties.F
    1274              :             CALL parallel_gemm('T', 'N', nmo_occ(ispin), nmo_virt(ispin), nao, 1.0_dp, &
    1275           10 :                                evects(ispin, istate), S_mos_virt(ispin), 0.0_dp, eigvec_X_ia_n(ispin))
    1276              : 
    1277              :             CALL get_exciton_descriptors(exc_descr, eigvec_X_ia_n(ispin), &
    1278              :                                          fm_multipole_ij, fm_multipole_ab, &
    1279              :                                          fm_multipole_ai, &
    1280           20 :                                          istate, nmo_occ(ispin), nmo_virt(ispin))
    1281              : 
    1282              :          END DO
    1283              :       END DO
    1284            2 :       CALL cp_fm_release(eigvec_X_ia_n)
    1285            2 :       CALL cp_fm_release(S_mos_virt)
    1286            2 :       CALL cp_fm_release(fm_multipole_ai)
    1287            2 :       CALL cp_fm_release(fm_multipole_ij)
    1288            2 :       CALL cp_fm_release(fm_multipole_ab)
    1289              : 
    1290              :       ! Actual printing
    1291            2 :       print_checkvalue = .TRUE.
    1292            2 :       prefix_output = ' '
    1293              :       CALL print_exciton_descriptors(exc_descr, ref_point_multipole, log_unit, &
    1294              :                                      nstates, print_checkvalue, do_directional_exciton_descriptors, &
    1295            2 :                                      prefix_output, qs_env)
    1296              : 
    1297            2 :       DEALLOCATE (ref_point_multipole)
    1298            2 :       DEALLOCATE (exc_descr)
    1299              : 
    1300            2 :       CALL timestop(handle)
    1301              : 
    1302            6 :    END SUBROUTINE tddfpt_print_exciton_descriptors
    1303              : 
    1304              : ! **************************************************************************************************
    1305              : !> \brief ...
    1306              : !> \param vin ...
    1307              : !> \param vout ...
    1308              : !> \param mos_occ ...
    1309              : !> \param matrix_s ...
    1310              : ! **************************************************************************************************
    1311            0 :    SUBROUTINE project_vector(vin, vout, mos_occ, matrix_s)
    1312              :       TYPE(dbcsr_type)                                   :: vin, vout
    1313              :       TYPE(cp_fm_type), INTENT(IN)                       :: mos_occ
    1314              :       TYPE(dbcsr_type), POINTER                          :: matrix_s
    1315              : 
    1316              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'project_vector'
    1317              : 
    1318              :       INTEGER                                            :: handle, nao, nmo
    1319              :       REAL(KIND=dp)                                      :: norm(1)
    1320              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_vec_struct
    1321              :       TYPE(cp_fm_type)                                   :: csvec, svec, vec
    1322              : 
    1323            0 :       CALL timeset(routineN, handle)
    1324              : 
    1325            0 :       CALL cp_fm_get_info(mos_occ, matrix_struct=fm_struct, nrow_global=nao, ncol_global=nmo)
    1326              :       CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
    1327            0 :                                nrow_global=nao, ncol_global=1)
    1328            0 :       CALL cp_fm_create(vec, fm_vec_struct)
    1329            0 :       CALL cp_fm_create(svec, fm_vec_struct)
    1330            0 :       CALL cp_fm_struct_release(fm_vec_struct)
    1331              :       CALL cp_fm_struct_create(fmstruct=fm_vec_struct, template_fmstruct=fm_struct, &
    1332            0 :                                nrow_global=nmo, ncol_global=1)
    1333            0 :       CALL cp_fm_create(csvec, fm_vec_struct)
    1334            0 :       CALL cp_fm_struct_release(fm_vec_struct)
    1335              : 
    1336            0 :       CALL copy_dbcsr_to_fm(vin, vec)
    1337            0 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, vec, svec, ncol=1, alpha=1.0_dp, beta=0.0_dp)
    1338            0 :       CALL parallel_gemm('T', 'N', nmo, 1, nao, 1.0_dp, mos_occ, svec, 0.0_dp, csvec)
    1339            0 :       CALL parallel_gemm('N', 'N', nao, 1, nmo, -1.0_dp, mos_occ, csvec, 1.0_dp, vec)
    1340            0 :       CALL cp_fm_vectorsnorm(vec, norm)
    1341            0 :       CPASSERT(norm(1) > 1.e-14_dp)
    1342            0 :       norm(1) = SQRT(1._dp/norm(1))
    1343            0 :       CALL cp_fm_scale(norm(1), vec)
    1344            0 :       CALL copy_fm_to_dbcsr(vec, vout, keep_sparsity=.FALSE.)
    1345              : 
    1346            0 :       CALL cp_fm_release(csvec)
    1347            0 :       CALL cp_fm_release(svec)
    1348            0 :       CALL cp_fm_release(vec)
    1349              : 
    1350            0 :       CALL timestop(handle)
    1351              : 
    1352            0 :    END SUBROUTINE project_vector
    1353              : 
    1354              : ! **************************************************************************************************
    1355              : !> \brief ...
    1356              : !> \param va ...
    1357              : !> \param vb ...
    1358              : !> \param res ...
    1359              : ! **************************************************************************************************
    1360            0 :    SUBROUTINE vec_product(va, vb, res)
    1361              :       TYPE(dbcsr_type)                                   :: va, vb
    1362              :       REAL(KIND=dp), INTENT(OUT)                         :: res
    1363              : 
    1364              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'vec_product'
    1365              : 
    1366              :       INTEGER                                            :: handle, icol, irow
    1367              :       LOGICAL                                            :: found
    1368            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vba, vbb
    1369              :       TYPE(dbcsr_iterator_type)                          :: iter
    1370              :       TYPE(mp_comm_type)                                 :: group
    1371              : 
    1372            0 :       CALL timeset(routineN, handle)
    1373              : 
    1374            0 :       res = 0.0_dp
    1375              : 
    1376            0 :       CALL dbcsr_get_info(va, group=group)
    1377            0 :       CALL dbcsr_iterator_start(iter, va)
    1378            0 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
    1379            0 :          CALL dbcsr_iterator_next_block(iter, irow, icol, vba)
    1380            0 :          CALL dbcsr_get_block_p(vb, row=irow, col=icol, block=vbb, found=found)
    1381            0 :          res = res + SUM(vba*vbb)
    1382            0 :          CPASSERT(found)
    1383              :       END DO
    1384            0 :       CALL dbcsr_iterator_stop(iter)
    1385            0 :       CALL group%sum(res)
    1386              : 
    1387            0 :       CALL timestop(handle)
    1388              : 
    1389            0 :    END SUBROUTINE vec_product
    1390              : 
    1391              : ! **************************************************************************************************
    1392              : !> \brief ...
    1393              : !> \param qs_env ...
    1394              : !> \param mos ...
    1395              : !> \param istate ...
    1396              : !> \param stride ...
    1397              : !> \param append_cube ...
    1398              : !> \param print_section ...
    1399              : ! **************************************************************************************************
    1400           16 :    SUBROUTINE print_nto_cubes(qs_env, mos, istate, stride, append_cube, print_section)
    1401              : 
    1402              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1403              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1404              :       INTEGER, INTENT(IN)                                :: istate
    1405              :       INTEGER, DIMENSION(:), POINTER                     :: stride
    1406              :       LOGICAL, INTENT(IN)                                :: append_cube
    1407              :       TYPE(section_vals_type), POINTER                   :: print_section
    1408              : 
    1409              :       CHARACTER(LEN=default_path_length)                 :: filename, my_pos_cube, title
    1410              :       INTEGER                                            :: i, iset, nmo, unit_nr
    1411              :       LOGICAL                                            :: mpi_io
    1412           16 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1413              :       TYPE(cell_type), POINTER                           :: cell
    1414              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1415              :       TYPE(cp_logger_type), POINTER                      :: logger
    1416              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1417              :       TYPE(particle_list_type), POINTER                  :: particles
    1418           16 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1419              :       TYPE(pw_c1d_gs_type)                               :: wf_g
    1420              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1421           16 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
    1422              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1423              :       TYPE(pw_r3d_rs_type)                               :: wf_r
    1424           16 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1425              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    1426              : 
    1427           32 :       logger => cp_get_default_logger()
    1428              : 
    1429           16 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
    1430           16 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
    1431           16 :       CALL auxbas_pw_pool%create_pw(wf_r)
    1432           16 :       CALL auxbas_pw_pool%create_pw(wf_g)
    1433              : 
    1434           16 :       CALL get_qs_env(qs_env, subsys=subsys)
    1435           16 :       CALL qs_subsys_get(subsys, particles=particles)
    1436              : 
    1437           16 :       my_pos_cube = "REWIND"
    1438           16 :       IF (append_cube) THEN
    1439            0 :          my_pos_cube = "APPEND"
    1440              :       END IF
    1441              : 
    1442              :       CALL get_qs_env(qs_env=qs_env, &
    1443              :                       atomic_kind_set=atomic_kind_set, &
    1444              :                       qs_kind_set=qs_kind_set, &
    1445              :                       cell=cell, &
    1446           16 :                       particle_set=particle_set)
    1447              : 
    1448           48 :       DO iset = 1, 2
    1449           32 :          CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
    1450           92 :          DO i = 1, nmo
    1451              :             CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
    1452           44 :                                         cell, dft_control, particle_set, pw_env)
    1453           44 :             IF (iset == 1) THEN
    1454           22 :                WRITE (filename, '(a4,I3.3,I2.2,a11)') "NTO_STATE", istate, i, "_Hole_State"
    1455           22 :             ELSEIF (iset == 2) THEN
    1456           22 :                WRITE (filename, '(a4,I3.3,I2.2,a15)') "NTO_STATE", istate, i, "_Particle_State"
    1457              :             END IF
    1458           44 :             mpi_io = .TRUE.
    1459              :             unit_nr = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
    1460              :                                            middle_name=TRIM(filename), file_position=my_pos_cube, &
    1461           44 :                                            log_filename=.FALSE., ignore_should_output=.TRUE., mpi_io=mpi_io)
    1462           44 :             IF (iset == 1) THEN
    1463           22 :                WRITE (title, *) "Natural Transition Orbital Hole State", i
    1464           22 :             ELSEIF (iset == 2) THEN
    1465           22 :                WRITE (title, *) "Natural Transition Orbital Particle State", i
    1466              :             END IF
    1467           44 :             CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=stride, mpi_io=mpi_io)
    1468              :             CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', &
    1469           76 :                                               ignore_should_output=.TRUE., mpi_io=mpi_io)
    1470              :          END DO
    1471              :       END DO
    1472              : 
    1473           16 :       CALL auxbas_pw_pool%give_back_pw(wf_g)
    1474           16 :       CALL auxbas_pw_pool%give_back_pw(wf_r)
    1475              : 
    1476           16 :    END SUBROUTINE print_nto_cubes
    1477              : 
    1478              : END MODULE qs_tddfpt2_properties
        

Generated by: LCOV version 2.0-1