61#include "./base/base_uses.f90"
67 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'bse_full_diag'
94 SUBROUTINE create_a(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
95 fm_A, Eigenval, unit_nr, &
96 homo, virtual, dimen_RI, mp2_env, &
99 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_s_ia_bse, fm_mat_s_bar_ij_bse, &
102 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenval
103 INTEGER,
INTENT(IN) :: unit_nr
104 INTEGER,
DIMENSION(:),
INTENT(IN) :: homo, virtual
105 INTEGER,
INTENT(IN) :: dimen_ri
106 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
110 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_A'
112 INTEGER :: a_virt_row, handle, i_occ_row, i_row_global, ii, isp, j_col_global, jj, k_isp, &
113 k_ov, n_ov_joint, ncol_local_a, nrow_local_a, nspins, sizeeigen
114 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: eig_offsets, n_ov, offsets
115 INTEGER,
DIMENSION(4) :: reordering
116 INTEGER,
DIMENSION(:),
POINTER :: col_indices_a, row_indices_a
117 REAL(kind=
dp) :: alpha, alpha_screening, eigen_diff
121 TYPE(
cp_fm_type) :: fm_a_copy, fm_s_joint, fm_w
126 CALL timeset(routinen, handle)
129 ALLOCATE (n_ov(nspins), offsets(nspins), eig_offsets(nspins))
134 eig_offsets(isp) = eig_offsets(isp - 1) + homo(isp - 1) + virtual(isp - 1)
137 NULLIFY (dft_control, tddfpt_control)
138 CALL get_qs_env(qs_env, dft_control=dft_control)
139 tddfpt_control => dft_control%tddfpt2_control
141 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print)
THEN
142 WRITE (unit_nr,
'(T2,A10,T13,A10)')
'BSE|DEBUG|',
'Creating A'
146 SELECT CASE (mp2_env%bse%bse_spin_config)
154 CALL cp_warn(__location__, &
155 "BSE: SPIN_CONFIG ignored for open-shell reference; using alpha=1.")
160 alpha_screening = mp2_env%bse%screening_factor
162 alpha_screening = 1.0_dp
175 nrow_global=n_ov_joint, ncol_global=n_ov_joint, &
176 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
180 IF (tddfpt_control%do_bse_w_only .AND. nspins == 1)
THEN
181 CALL cp_fm_create(fm_a_copy, fm_struct_a, name=
"fm_A_iajb")
187 IF ((.NOT. tddfpt_control%do_bse) .AND. (.NOT. tddfpt_control%do_bse_w_only))
THEN
191 context=fm_mat_s_ia_bse(1)%matrix_struct%context, &
192 nrow_global=dimen_ri, ncol_global=n_ov_joint, &
193 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
194 CALL cp_fm_create(fm_s_joint, fm_struct_s_joint, name=
"fm_S_ia_joint")
197 CALL parallel_gemm(transa=
"T", transb=
"N", m=n_ov_joint, n=n_ov_joint, k=dimen_ri, &
198 alpha=alpha, matrix_a=fm_s_joint, matrix_b=fm_s_joint, beta=0.0_dp, &
203 CALL parallel_gemm(transa=
"T", transb=
"N", m=homo(1)*virtual(1), n=homo(1)*virtual(1), &
204 k=dimen_ri, alpha=alpha, &
205 matrix_a=fm_mat_s_ia_bse(1), matrix_b=fm_mat_s_ia_bse(1), &
206 beta=0.0_dp, matrix_c=fm_a)
210 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print)
THEN
211 WRITE (unit_nr,
'(T2,A10,T13,A16)')
'BSE|DEBUG|',
'Allocated A_iajb'
219 CALL cp_fm_struct_create(fm_struct_w, context=fm_mat_s_ab_bse(isp)%matrix_struct%context, &
220 nrow_global=homo(isp)**2, ncol_global=virtual(isp)**2, &
221 para_env=fm_mat_s_ab_bse(isp)%matrix_struct%para_env)
225 CALL parallel_gemm(transa=
"T", transb=
"N", m=homo(isp)**2, n=virtual(isp)**2, &
226 k=dimen_ri, alpha=alpha_screening, &
227 matrix_a=fm_mat_s_bar_ij_bse(isp), matrix_b=fm_mat_s_ab_bse(isp), &
228 beta=0.0_dp, matrix_c=fm_w)
229 reordering = [1, 3, 2, 4]
231 virtual(isp), virtual(isp), unit_nr, reordering, mp2_env, &
232 row_offset=offsets(isp), col_offset=offsets(isp))
233 IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only)
THEN
235 virtual(1), virtual(1), unit_nr, reordering, mp2_env)
238 IF (nspins == 1)
THEN
239 IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
240 tddfpt_control%do_bse_gw_only)
THEN
243 IF (.NOT. tddfpt_control%do_bse_gw_only)
THEN
244 ALLOCATE (ex_env%bse_w_matrix_MO(1, 1))
245 ALLOCATE (ex_env%bse_a_matrix_MO(1, 1))
246 CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_w)
247 CALL cp_fm_create(ex_env%bse_a_matrix_MO(1, 1), fm_struct_a)
248 CALL cp_fm_to_fm(fm_w, ex_env%bse_w_matrix_MO(1, 1))
249 IF (tddfpt_control%do_bse_w_only)
THEN
250 CALL cp_fm_to_fm(fm_a_copy, ex_env%bse_a_matrix_MO(1, 1))
252 CALL cp_fm_to_fm(fm_a, ex_env%bse_a_matrix_MO(1, 1))
262 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print)
THEN
263 WRITE (unit_nr,
'(T2,A10,T13,A16)')
'BSE|DEBUG|',
'Allocated W_ijab'
265 IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only)
CALL cp_fm_release(fm_a_copy)
268 CALL cp_fm_get_info(matrix=fm_a, nrow_local=nrow_local_a, ncol_local=ncol_local_a, &
269 row_indices=row_indices_a, col_indices=col_indices_a)
272 IF (.NOT. tddfpt_control%do_bse)
THEN
273 DO ii = 1, nrow_local_a
274 i_row_global = row_indices_a(ii)
275 DO jj = 1, ncol_local_a
276 j_col_global = col_indices_a(jj)
277 IF (i_row_global == j_col_global)
THEN
280 DO k_isp = 1, nspins - 1
281 IF (i_row_global <= offsets(k_isp) + n_ov(k_isp))
THEN
286 k_ov = i_row_global - offsets(isp)
287 i_occ_row = (k_ov - 1)/virtual(isp) + 1
288 a_virt_row = mod(k_ov - 1, virtual(isp)) + 1
289 eigen_diff = eigenval(eig_offsets(isp) + a_virt_row + homo(isp)) - &
290 eigenval(eig_offsets(isp) + i_occ_row)
291 fm_a%local_data(ii, jj) = fm_a%local_data(ii, jj) + eigen_diff
298 IF (nspins == 1)
THEN
299 IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
300 tddfpt_control%do_bse_gw_only)
THEN
301 sizeeigen =
SIZE(eigenval)
302 ALLOCATE (ex_env%gw_eigen(sizeeigen))
303 ex_env%gw_eigen(:) = eigenval(:)
308 DEALLOCATE (n_ov, offsets, eig_offsets)
312 CALL timestop(handle)
331 SUBROUTINE create_b(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &
332 homo, virtual, dimen_RI, unit_nr, mp2_env)
334 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mat_s_ia_bse, fm_mat_s_bar_ia_bse
336 INTEGER,
DIMENSION(:),
INTENT(IN) :: homo, virtual
337 INTEGER,
INTENT(IN) :: dimen_ri, unit_nr
338 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
340 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_B'
342 INTEGER :: handle, isp, n_ov_joint, nspins
343 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: n_ov, offsets
344 INTEGER,
DIMENSION(4) :: reordering
345 REAL(kind=
dp) :: alpha, alpha_screening
350 CALL timeset(routinen, handle)
353 ALLOCATE (n_ov(nspins), offsets(nspins))
356 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print)
THEN
357 WRITE (unit_nr,
'(T2,A10,T13,A10)')
'BSE|DEBUG|',
'Creating B'
362 SELECT CASE (mp2_env%bse%bse_spin_config)
368 IF (nspins > 1) alpha = 1.0_dp
371 alpha_screening = mp2_env%bse%screening_factor
373 alpha_screening = 1.0_dp
377 NULLIFY (fm_struct_b)
379 nrow_global=n_ov_joint, ncol_global=n_ov_joint, &
380 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
386 NULLIFY (fm_struct_s_joint)
388 context=fm_mat_s_ia_bse(1)%matrix_struct%context, &
389 nrow_global=dimen_ri, ncol_global=n_ov_joint, &
390 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
391 CALL cp_fm_create(fm_s_joint, fm_struct_s_joint, name=
"fm_S_ia_joint")
394 CALL parallel_gemm(transa=
"T", transb=
"N", m=n_ov_joint, n=n_ov_joint, k=dimen_ri, &
395 alpha=alpha, matrix_a=fm_s_joint, matrix_b=fm_s_joint, beta=0.0_dp, &
400 CALL parallel_gemm(transa=
"T", transb=
"N", m=homo(1)*virtual(1), n=homo(1)*virtual(1), &
401 k=dimen_ri, alpha=alpha, &
402 matrix_a=fm_mat_s_ia_bse(1), matrix_b=fm_mat_s_ia_bse(1), &
403 beta=0.0_dp, matrix_c=fm_b)
406 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print)
THEN
407 WRITE (unit_nr,
'(T2,A10,T13,A16)')
'BSE|DEBUG|',
'Allocated B_iajb'
414 NULLIFY (fm_struct_w)
416 context=fm_mat_s_ia_bse(isp)%matrix_struct%context, &
417 nrow_global=homo(isp)*virtual(isp), &
418 ncol_global=homo(isp)*virtual(isp), &
419 para_env=fm_mat_s_ia_bse(isp)%matrix_struct%para_env)
422 CALL parallel_gemm(transa=
"T", transb=
"N", m=homo(isp)*virtual(isp), &
423 n=homo(isp)*virtual(isp), k=dimen_ri, alpha=alpha_screening, &
424 matrix_a=fm_mat_s_bar_ia_bse(isp), matrix_b=fm_mat_s_ia_bse(isp), &
425 beta=0.0_dp, matrix_c=fm_w)
426 reordering = [1, 4, 3, 2]
428 virtual(isp), virtual(isp), unit_nr, reordering, mp2_env, &
429 row_offset=offsets(isp), col_offset=offsets(isp))
436 DEALLOCATE (n_ov, offsets)
438 CALL timestop(handle)
457 fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
458 unit_nr, mp2_env, diag_est)
461 TYPE(
cp_fm_type),
INTENT(INOUT) :: fm_c, fm_sqrt_a_minus_b, &
462 fm_inv_sqrt_a_minus_b
463 INTEGER,
INTENT(IN) :: unit_nr
464 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
465 REAL(kind=
dp),
INTENT(IN) :: diag_est
467 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_hermitian_form_of_ABBA'
469 INTEGER :: dim_mat, handle, n_dependent
470 REAL(kind=
dp),
DIMENSION(2) :: eigvals_ab_diff
471 TYPE(
cp_fm_type) :: fm_a_minus_b, fm_a_plus_b, fm_dummy, &
474 CALL timeset(routinen, handle)
476 IF (unit_nr > 0)
THEN
477 WRITE (unit_nr,
'(T2,A4,T7,A25,A39,ES6.0,A3)')
'BSE|',
'Diagonalizing aux. matrix', &
478 ' with size of A. This will take around ', diag_est,
" s."
495 CALL cp_fm_create(fm_sqrt_a_minus_b, fm_a%matrix_struct)
497 CALL cp_fm_create(fm_inv_sqrt_a_minus_b, fm_a%matrix_struct)
502 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print)
THEN
503 WRITE (unit_nr,
'(T2,A10,T13,A19)')
'BSE|DEBUG|',
'Created work arrays'
511 CALL cp_fm_to_fm(fm_a_minus_b, fm_inv_sqrt_a_minus_b)
519 CALL cp_fm_power(fm_inv_sqrt_a_minus_b, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_ab_diff)
523 IF (eigvals_ab_diff(1) < 0)
THEN
524 CALL cp_abort(__location__, &
525 "Matrix (A-B) is not positive definite. "// &
526 "Hermitian diagonalization of full ABBA matrix is ill-defined.")
533 CALL parallel_gemm(
"N",
"N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_inv_sqrt_a_minus_b, fm_a_minus_b, 0.0_dp, &
537 CALL parallel_gemm(
"N",
"N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_sqrt_a_minus_b, fm_a_plus_b, 0.0_dp, &
548 CALL parallel_gemm(
"N",
"N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_work_product, fm_sqrt_a_minus_b, 0.0_dp, &
552 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print)
THEN
553 WRITE (unit_nr,
'(T2,A10,T13,A36)')
'BSE|DEBUG|',
'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
556 CALL timestop(handle)
576 fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
577 unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
580 INTEGER,
DIMENSION(:),
INTENT(IN) :: homo, virtual, homo_irred
581 TYPE(
cp_fm_type),
INTENT(INOUT) :: fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b
582 INTEGER,
INTENT(IN) :: unit_nr
583 REAL(kind=
dp),
INTENT(IN) :: diag_est
584 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
586 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: mo_coeff
588 CHARACTER(LEN=*),
PARAMETER :: routinen =
'diagonalize_C'
590 INTEGER :: diag_info, handle, n_ov_joint, nspins
591 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: exc_ens
592 TYPE(
cp_fm_type) :: fm_eigvec_x, fm_eigvec_y, fm_eigvec_z, &
593 fm_mat_eigvec_transform_diff, &
594 fm_mat_eigvec_transform_sum
596 CALL timeset(routinen, handle)
599 n_ov_joint = sum(homo*virtual)
601 IF (unit_nr > 0)
THEN
602 WRITE (unit_nr,
'(T2,A4,T7,A17,A22,ES6.0,A3)')
'BSE|',
'Diagonalizing C. ', &
603 'This will take around ', diag_est,
' s.'
610 ALLOCATE (exc_ens(n_ov_joint))
614 IF (diag_info /= 0)
THEN
615 CALL cp_abort(__location__, &
616 "Diagonalization of C=(A-B)^0.5 (A+B) (A-B)^0.5 failed in BSE")
622 IF (exc_ens(1) < 0)
THEN
623 IF (unit_nr > 0)
THEN
624 CALL cp_abort(__location__, &
625 "Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
626 "(A+B) is not positive definite.")
629 exc_ens = sqrt(exc_ens)
640 CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_c%matrix_struct)
642 CALL parallel_gemm(transa=
"N", transb=
"N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
643 matrix_a=fm_sqrt_a_minus_b, matrix_b=fm_eigvec_z, beta=0.0_dp, &
644 matrix_c=fm_mat_eigvec_transform_sum)
650 CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_c%matrix_struct)
652 CALL parallel_gemm(transa=
"N", transb=
"N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
653 matrix_a=fm_inv_sqrt_a_minus_b, matrix_b=fm_eigvec_z, beta=0.0_dp, &
654 matrix_c=fm_mat_eigvec_transform_diff)
664 CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_x)
670 CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_y)
677 IF (nspins == 1)
THEN
678 CALL postprocess_bse(exc_ens, fm_eigvec_x, mp2_env, qs_env, mo_coeff, &
679 homo(1), virtual(1), homo_irred(1), unit_nr, &
680 .false., fm_eigvec_y)
683 CALL bse_open_shell_optical(exc_ens, fm_eigvec_x, homo, virtual, homo_irred, &
684 .false., qs_env, mo_coeff, mp2_env, unit_nr, fm_eigvec_y)
691 CALL timestop(handle)
708 unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
711 INTEGER,
DIMENSION(:),
INTENT(IN) :: homo, virtual, homo_irred
712 INTEGER,
INTENT(IN) :: unit_nr
713 REAL(kind=
dp),
INTENT(IN) :: diag_est
714 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
716 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: mo_coeff
718 CHARACTER(LEN=*),
PARAMETER :: routinen =
'diagonalize_A'
720 INTEGER :: diag_info, handle, n_ov_joint, nspins
721 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: exc_ens
724 CALL timeset(routinen, handle)
727 n_ov_joint = sum(homo*virtual)
729 IF (unit_nr > 0)
THEN
730 WRITE (unit_nr,
'(T2,A4,T7,A17,A22,ES6.0,A3)')
'BSE|',
'Diagonalizing A. ', &
731 'This will take around ', diag_est,
' s.'
736 ALLOCATE (exc_ens(n_ov_joint))
740 IF (diag_info /= 0)
THEN
741 CALL cp_abort(__location__, &
742 "Diagonalization of A failed in TDA-BSE")
745 IF (nspins == 1)
THEN
746 CALL postprocess_bse(exc_ens, fm_eigvec, mp2_env, qs_env, mo_coeff, &
747 homo(1), virtual(1), homo_irred(1), unit_nr, .true.)
749 CALL bse_open_shell_optical(exc_ens, fm_eigvec, homo, virtual, homo_irred, &
750 .true., qs_env, mo_coeff, mp2_env, unit_nr)
756 CALL timestop(handle)
776 SUBROUTINE bse_open_shell_optical(Exc_ens, fm_eigvec_X, homo, virtual, homo_irred, &
777 flag_tda, qs_env, mo_coeff, mp2_env, unit_nr, fm_eigvec_Y)
779 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: exc_ens
781 INTEGER,
DIMENSION(:),
INTENT(IN) :: homo, virtual, homo_irred
782 LOGICAL,
INTENT(IN) :: flag_tda
784 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: mo_coeff
785 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
786 INTEGER,
INTENT(IN) :: unit_nr
787 TYPE(
cp_fm_type),
INTENT(IN),
OPTIONAL :: fm_eigvec_y
789 CHARACTER(LEN=*),
PARAMETER :: routinen =
'bse_open_shell_optical'
791 CHARACTER(LEN=10) :: info_approximation, multiplet
792 INTEGER :: handle, idir, isp, jdir, n, n_ov_joint, &
794 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: n_ov_sp, offsets_sp
795 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: oscill_str_joint, ref_pt
796 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: pol_res_joint, trans_mom_joint
799 TYPE(
cp_fm_type) :: fm_dip_reord_sp, fm_eigvec_sp, &
801 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: fm_dip_ab_sp, fm_dip_ai_sp, fm_dip_ij_sp
802 TYPE(
cp_fm_type),
DIMENSION(3) :: fm_trans_mom_joint
804 CALL timeset(routinen, handle)
807 ALLOCATE (n_ov_sp(nspins), offsets_sp(nspins))
814 info_approximation =
" -TDA- "
816 info_approximation =
"-ABBA-"
819 IF (unit_nr > 0)
THEN
820 WRITE (unit_nr,
'(T2,A4,T7,A43)')
'BSE|',
'Joint open-shell BSE excitation energies:'
823 info_approximation, mp2_env, unit_nr)
827 info_approximation, mp2_env, unit_nr, fm_eigvec_y)
830 CALL cp_fm_create(fm_trans_coeff, fm_eigvec_x%matrix_struct)
832 IF (
PRESENT(fm_eigvec_y))
CALL cp_fm_scale_and_add(1.0_dp, fm_trans_coeff, 1.0_dp, fm_eigvec_y)
836 ALLOCATE (fm_dip_ai_sp(3), fm_dip_ij_sp(3), fm_dip_ab_sp(3), ref_pt(3))
837 ALLOCATE (oscill_str_joint(n_ov_joint), trans_mom_joint(3, 1, n_ov_joint))
838 ALLOCATE (pol_res_joint(3, 3, n_ov_joint))
839 trans_mom_joint(:, :, :) = 0.0_dp
840 NULLIFY (fm_struct_dip_reord, fm_struct_sp, fm_struct_tmom)
842 fm_trans_coeff%matrix_struct%context, 1, n_ov_joint)
844 CALL cp_fm_create(fm_trans_mom_joint(idir), fm_struct_tmom)
849 qs_env, mo_coeff(isp:isp), ref_pt, 1, &
850 homo(isp), virtual(isp), fm_trans_coeff%matrix_struct%context, &
852 NULLIFY (fm_struct_sp, fm_struct_dip_reord)
854 fm_trans_coeff%matrix_struct%context, n_ov_sp(isp), n_ov_joint)
858 offsets_sp(isp) + 1, 1, 1, 1)
860 fm_trans_coeff%matrix_struct%context, 1, n_ov_sp(isp))
862 CALL cp_fm_create(fm_dip_reord_sp, fm_struct_dip_reord, name=
"bse_dip_reord")
865 1, 1, 1, virtual(isp), unit_nr, [2, 4, 3, 1], mp2_env)
866 CALL parallel_gemm(
'N',
'N', 1, n_ov_joint, n_ov_sp(isp), 1.0_dp, &
867 fm_dip_reord_sp, fm_eigvec_sp, 1.0_dp, fm_trans_mom_joint(idir))
875 NULLIFY (fm_struct_sp)
877 NULLIFY (fm_struct_dip_reord)
887 pol_res_joint(idir, jdir, n) = 2.0_dp*exc_ens(n)*trans_mom_joint(idir, 1, n) &
888 *trans_mom_joint(jdir, 1, n)
891 oscill_str_joint(n) = 2.0_dp/3.0_dp*exc_ens(n)*sum(abs(trans_mom_joint(:, 1, n))**2)
894 n_ov_joint, 1, sum(homo_irred), flag_tda, info_approximation, &
895 mp2_env, unit_nr, open_shell=.true.)
897 CALL cp_warn(__location__, &
898 "Open-shell (UKS) BSE: exciton descriptors and NTO analysis are not yet "// &
899 "implemented and have been skipped.")
901 DEALLOCATE (fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp, ref_pt)
902 DEALLOCATE (n_ov_sp, offsets_sp, oscill_str_joint, trans_mom_joint, pol_res_joint)
904 CALL timestop(handle)
906 END SUBROUTINE bse_open_shell_optical
922 SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
923 homo, virtual, homo_irred, unit_nr, &
924 flag_TDA, fm_eigvec_Y)
926 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: exc_ens
928 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
930 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: mo_coeff
931 INTEGER :: homo, virtual, homo_irred, unit_nr
932 LOGICAL,
OPTIONAL :: flag_tda
933 TYPE(
cp_fm_type),
INTENT(IN),
OPTIONAL :: fm_eigvec_y
935 CHARACTER(LEN=*),
PARAMETER :: routinen =
'postprocess_bse'
937 CHARACTER(LEN=10) :: info_approximation, multiplet
938 INTEGER :: handle, i_exc, idir, n_moments_di, &
940 REAL(kind=
dp) :: alpha
941 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: oscill_str, ref_point_multipole
942 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
944 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
945 fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
947 DIMENSION(:) :: exc_descr
949 CALL timeset(routinen, handle)
952 IF (mp2_env%bse%bse_spin_config == 0)
THEN
953 multiplet =
"Singlet"
956 multiplet =
"Triplet"
959 IF (.NOT.
PRESENT(flag_tda))
THEN
963 info_approximation =
" -TDA- "
965 info_approximation =
"-ABBA-"
972 ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
973 ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
974 ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
975 ALLOCATE (ref_point_multipole(3))
977 CALL get_multipoles_mo(fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
978 qs_env, mo_coeff, ref_point_multipole, 1, &
979 homo, virtual, fm_eigvec_x%matrix_struct%context)
981 IF (mp2_env%bse%num_print_exc_descr > 0)
THEN
983 ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
984 ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
985 ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
986 CALL get_multipoles_mo(fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
987 qs_env, mo_coeff, ref_point_multipole, 2, &
988 homo, virtual, fm_eigvec_x%matrix_struct%context)
990 ALLOCATE (exc_descr(mp2_env%bse%num_print_exc_descr))
991 DO i_exc = 1, mp2_env%bse%num_print_exc_descr
993 .false., unit_nr, mp2_env)
994 IF (.NOT. flag_tda)
THEN
996 .false., unit_nr, mp2_env)
999 fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1000 fm_quadpole_ai_trunc, &
1001 i_exc, homo, virtual, &
1005 fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1006 fm_quadpole_ai_trunc, &
1007 i_exc, homo, virtual)
1010 IF (.NOT. flag_tda)
THEN
1016 IF (mp2_env%bse%bse_spin_config == 0)
THEN
1018 trans_mom_bse, oscill_str, polarizability_residues, &
1019 mp2_env, homo, virtual, unit_nr, &
1025 multiplet, alpha, mp2_env, unit_nr)
1029 info_approximation, mp2_env, unit_nr)
1033 info_approximation, mp2_env, unit_nr, fm_eigvec_y)
1037 homo, virtual, homo_irred, flag_tda, &
1038 info_approximation, mp2_env, unit_nr)
1040 IF (mp2_env%bse%num_print_exc_descr > 0)
THEN
1042 mp2_env%bse%num_print_exc_descr, mp2_env%bse%bse_debug_print, &
1043 mp2_env%bse%print_directional_exc_descr, &
1048 IF (mp2_env%bse%do_nto_analysis)
THEN
1049 IF (unit_nr > 0)
THEN
1050 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
1051 WRITE (unit_nr,
'(T2,A4,T7,A47)') &
1052 'BSE|',
"Calculating Natural Transition Orbitals (NTOs)."
1053 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
1056 mo_coeff, homo, virtual, &
1057 info_approximation, &
1059 qs_env, unit_nr, mp2_env)
1062 DO idir = 1, n_moments_di
1067 IF (mp2_env%bse%num_print_exc_descr > 0)
THEN
1068 DO idir = 1, n_moments_quad
1073 DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
1074 DEALLOCATE (exc_descr)
1076 DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
1077 DEALLOCATE (ref_point_multipole)
1078 IF (mp2_env%bse%bse_spin_config == 0)
THEN
1079 DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
1081 IF (unit_nr > 0)
THEN
1082 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
1083 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
1086 CALL timestop(handle)
1088 END SUBROUTINE postprocess_bse
Routines for the full diagonalization of GW + Bethe-Salpeter for computing electronic excitations.
subroutine, public diagonalize_a(fm_a, homo, virtual, homo_irred, unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
Solving hermitian eigenvalue equation A X^n = Ω^n X^n.
subroutine, public diagonalize_c(fm_c, homo, virtual, homo_irred, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n . Here, the eigenvectors Z^n relate to X^n via Eq....
subroutine, public create_hermitian_form_of_abba(fm_a, fm_b, fm_c, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, mp2_env, diag_est)
Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem (cf....
subroutine, public create_b(fm_mat_s_ia_bse, fm_mat_s_bar_ia_bse, fm_b, homo, virtual, dimen_ri, unit_nr, mp2_env)
Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W) B_ia,jb = α * v_ia,...
subroutine, public create_a(fm_mat_s_ia_bse, fm_mat_s_bar_ij_bse, fm_mat_s_ab_bse, fm_a, eigenval, unit_nr, homo, virtual, dimen_ri, mp2_env, para_env, qs_env)
Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W) A_ia,...
Routines for printing information in context of the BSE calculation.
subroutine, public print_output_header(homo, virtual, homo_irred, flag_tda, multiplet, alpha, mp2_env, unit_nr)
...
subroutine, public print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, num_print_exc_descr, print_checkvalue, print_directional_exc_descr, prefix_output, qs_env)
Prints exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018).
subroutine, public print_optical_properties(exc_ens, oscill_str, trans_mom_bse, polarizability_residues, homo, virtual, homo_irred, flag_tda, info_approximation, mp2_env, unit_nr, open_shell)
...
subroutine, public print_transition_amplitudes(fm_eigvec_x, homo, virtual, homo_irred, info_approximation, mp2_env, unit_nr, fm_eigvec_y)
...
subroutine, public print_excitation_energies(exc_ens, homo, virtual, flag_tda, multiplet, info_approximation, mp2_env, unit_nr)
...
Routines for computing excitonic properties, e.g. exciton diameter, from the BSE.
subroutine, public get_oscillator_strengths(fm_eigvec_x, exc_ens, fm_dipole_ai_trunc, trans_mom_bse, oscill_str, polarizability_residues, mp2_env, homo_red, virtual_red, unit_nr, fm_eigvec_y)
Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n ) and oscil...
subroutine, public get_exciton_descriptors(exc_descr, fm_x_ia, fm_multipole_ij_trunc, fm_multipole_ab_trunc, fm_multipole_ai_trunc, i_exc, homo, virtual, fm_y_ia)
...
subroutine, public calculate_ntos(fm_x, fm_y, mo_coeff, homo, virtual, info_approximation, oscill_str, qs_env, unit_nr, mp2_env)
...
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
subroutine, public assemble_joint_ov_slab(fm_s_ia, offsets, n_ov, dimen_ri, fm_s_joint)
Column-concatenate per-spin ia-slabs into the joint dimen_RI x n_ov_joint slab. Sigma-block of spin i...
subroutine, public fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env, row_offset, col_offset)
Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,...
subroutine, public get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, qs_env, mo_coeff, rpoint, n_moments, homo_red, virtual_red, context_bse, ispin)
...
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, unit_nr, mp2_env)
...
subroutine, public get_bse_spin_block_layout(homo_red, virt_red, n_ov, offsets, n_ov_joint)
Spin-block layout for the open-shell (joint) BSE matrix: per-spin OV-pair counts and the block offset...
subroutine, public comp_eigvec_coeff_bse(fm_work, eig_vals, beta, gamma, do_transpose)
Routine for computing the coefficients of the eigenvectors of the BSE matrix from a multiplication wi...
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Types for excited states potential energies.
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Defines the basic variable types.
integer, parameter, public dp
Interface to the message passing library MPI.
Types needed for MP2 calculations.
basic linear algebra operations for full matrixes
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
Contains information on the excited states energy.
stores all the informations relevant to an mpi environment