52#include "./base/base_uses.f90"
58 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'bse_properties'
68 REAL(kind=
dp),
DIMENSION(3) :: r_e = 0.0_dp, &
75 sigma_e_dir = 0.0_dp, &
76 sigma_h_dir = 0.0_dp, &
78 REAL(kind=
dp),
DIMENSION(3, 3) :: r_e_h = 0.0_dp, &
80 corr_e_h_matrix = 0.0_dp
81 REAL(kind=
dp) :: sigma_e = 0.0_dp, &
83 cov_e_h_sum = 0.0_dp, &
85 diff_r_abs = 0.0_dp, &
86 diff_r_sqr = 0.0_dp, &
88 LOGICAL :: flag_tda = .false.
112 trans_mom_bse, oscill_str, polarizability_residues, &
113 mp2_env, homo_red, virtual_red, unit_nr, &
117 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
118 INTENT(IN) :: exc_ens
119 TYPE(
cp_fm_type),
DIMENSION(3) :: fm_dipole_ai_trunc
120 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
121 INTENT(OUT) :: trans_mom_bse
122 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
123 INTENT(OUT) :: oscill_str
124 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
125 INTENT(OUT) :: polarizability_residues
126 TYPE(
mp2_type),
INTENT(IN) :: mp2_env
127 INTEGER,
INTENT(IN) :: homo_red, virtual_red, unit_nr
128 TYPE(
cp_fm_type),
INTENT(IN),
OPTIONAL :: fm_eigvec_y
130 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_oscillator_strengths'
132 INTEGER :: handle, idir, jdir, n
134 fm_struct_trans_mom_bse
136 TYPE(
cp_fm_type),
DIMENSION(3) :: fm_dipole_mo_trunc_reordered, &
137 fm_dipole_per_dir, fm_trans_mom_bse
139 CALL timeset(routinen, handle)
141 CALL cp_fm_struct_create(fm_struct_dipole_mo_trunc_reordered, fm_eigvec_x%matrix_struct%para_env, &
142 fm_eigvec_x%matrix_struct%context, 1, homo_red*virtual_red)
144 fm_eigvec_x%matrix_struct%context, 1, homo_red*virtual_red)
151 CALL cp_fm_create(fm_dipole_mo_trunc_reordered(idir), matrix_struct=fm_struct_dipole_mo_trunc_reordered, &
152 name=
"dipoles_mo_reordered")
153 CALL cp_fm_set_all(fm_dipole_mo_trunc_reordered(idir), 0.0_dp)
154 CALL fm_general_add_bse(fm_dipole_mo_trunc_reordered(idir), fm_dipole_ai_trunc(idir), 1.0_dp, &
157 unit_nr, [2, 4, 3, 1], mp2_env)
162 CALL cp_fm_create(fm_trans_mom_bse(idir), matrix_struct=fm_struct_trans_mom_bse, &
163 name=
"excitonic_dipoles")
168 CALL cp_fm_create(fm_eigvec_xysum, matrix_struct=fm_eigvec_x%matrix_struct, name=
"excit_amplitude_sum")
171 IF (
PRESENT(fm_eigvec_y))
THEN
175 CALL parallel_gemm(
'N',
'N', 1, homo_red*virtual_red, homo_red*virtual_red, sqrt(2.0_dp), &
176 fm_dipole_mo_trunc_reordered(idir), fm_eigvec_xysum, 0.0_dp, fm_trans_mom_bse(idir))
180 ALLOCATE (oscill_str(homo_red*virtual_red))
183 ALLOCATE (trans_mom_bse(3, 1, homo_red*virtual_red))
184 ALLOCATE (polarizability_residues(3, 3, homo_red*virtual_red))
185 trans_mom_bse(:, :, :) = 0.0_dp
192 DO n = 1, homo_red*virtual_red
195 polarizability_residues(idir, jdir, n) = 2.0_dp*exc_ens(n)*trans_mom_bse(idir, 1, n)*trans_mom_bse(jdir, 1, n)
198 oscill_str(n) = 2.0_dp/3.0_dp*exc_ens(n)*sum(abs(trans_mom_bse(:, 1, n))**2)
210 CALL timestop(handle)
227 info_approximation, unit_nr, mp2_env)
229 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
230 INTENT(IN) :: oscill_str
231 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :), &
232 INTENT(IN) :: polarizability_residues
233 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
234 INTENT(IN) :: exc_ens
235 CHARACTER(LEN=10) :: info_approximation
236 INTEGER,
INTENT(IN) :: unit_nr
237 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
239 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_and_print_absorption_spectrum'
241 CHARACTER(LEN=10) :: eta_str, width_eta_format_str
242 CHARACTER(LEN=40) :: file_name_crosssection, &
244 INTEGER :: handle, i, idir, j, jdir, k, num_steps, &
245 unit_nr_file, width_eta
246 REAL(kind=
dp) :: eta, freq_end, freq_start, freq_step, &
248 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: abs_cross_section, abs_spectrum
249 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eta_list
252 CALL timeset(routinen, handle)
254 freq_step = mp2_env%bse%bse_spectrum_freq_step_size
255 freq_start = mp2_env%bse%bse_spectrum_freq_start
256 freq_end = mp2_env%bse%bse_spectrum_freq_end
257 eta_list => mp2_env%bse%bse_eta_spectrum_list
260 num_steps = nint((freq_end - freq_start)/freq_step) + 1
262 DO k = 1,
SIZE(eta_list)
266 width_eta = max(1, int(log10(eta)) + 1) + 4
267 WRITE (width_eta_format_str,
"(A2,I0,A3)")
'(F', width_eta,
'.3)'
268 WRITE (eta_str, width_eta_format_str) eta*
evolt
270 file_name_spectrum =
'BSE'//trim(adjustl(info_approximation))//
'eta='//trim(eta_str)//
'.spectrum'
271 file_name_crosssection =
'BSE'//trim(adjustl(info_approximation))//
'eta='//trim(eta_str)//
'.crosssection'
275 ALLOCATE (abs_spectrum(num_steps, 11))
276 abs_spectrum(:, :) = 0.0_dp
279 ALLOCATE (abs_cross_section(num_steps, 11))
280 abs_cross_section(:, :) = 0.0_dp
295 omega = freq_start + (i - 1)*freq_step
296 abs_spectrum(i, 1) = omega
297 DO j = 1,
SIZE(oscill_str)
298 abs_spectrum(i, 2) = abs_spectrum(i, 2) - oscill_str(j)* &
299 aimag(1/((omega + cmplx(0.0, eta, kind=
dp))**2 - exc_ens(j)**2))
304 abs_spectrum(i, 2 + (idir - 1)*3 + jdir) = abs_spectrum(i, 2 + (idir - 1)*3 + jdir) &
305 - polarizability_residues(idir, jdir, j)* &
306 aimag(1/((omega + cmplx(0.0, eta, kind=
dp))**2 - exc_ens(j)**2))
314 omega = abs_spectrum(i, 1)
315 abs_cross_section(i, 1) = omega
316 abs_cross_section(i, 2:) = 4.0_dp*
pi*abs_spectrum(i, 2:)*omega/
c_light_au
320 IF (mp2_env%bse%bse_debug_print)
THEN
321 IF (unit_nr > 0)
THEN
322 WRITE (unit_nr,
'(T2,A10,T13,A,T65,F16.4)')
'BSE|DEBUG|', &
323 'Averaged dynamical dipole polarizability at 8.2 eV:', &
325 WRITE (unit_nr,
'(T2,A10,T13,A,T65,F16.4)')
'BSE|DEBUG|', &
326 'Averaged photoabsorption cross section at 8.2 eV:', &
327 abs_cross_section(83, 2)
333 IF (logger%para_env%is_source())
THEN
339 IF (unit_nr_file > 0)
THEN
340 CALL open_file(file_name_crosssection, unit_number=unit_nr_file, &
341 file_status=
"UNKNOWN", file_action=
"WRITE")
342 WRITE (unit_nr_file,
'(A,A6)') σµµωπω
"# Photoabsorption cross section _{ '}() = -4/c * Im[ \sum_n "// &
343 Ωµµωη²Ω²
"[2 ^n d_^n d_'^n] / [(+i)- (^n)] ] from Bethe Salpeter equation for method ", &
344 trim(adjustl(info_approximation))
345 WRITE (unit_nr_file,
'(A20,1X,10(2X,A20,1X))')
"# Frequency (eV)", σω
"_{avg}()", σω
"_xx()", &
346 σω
"_xy()", σω
"_xz()", σω
"_yx()", σω
"_yy()", σω
"_yz()", σω
"_zx()", &
349 WRITE (unit_nr_file,
'(11(F20.8,1X))') abs_cross_section(i, 1)*
evolt, abs_cross_section(i, 2:11)
353 DEALLOCATE (abs_cross_section)
355 IF (unit_nr_file > 0)
THEN
356 CALL open_file(file_name_spectrum, unit_number=unit_nr_file, &
357 file_status=
"UNKNOWN", file_action=
"WRITE")
358 WRITE (unit_nr_file,
'(A,A6)') αµµω
"# Imaginary part of polarizability _{ '}() = -\sum_n "// &
359 Ωµµωη²Ω²
"[2 ^n d_^n d_'^n] / [(+i)- (^n)] from Bethe Salpeter equation for method ", &
360 trim(adjustl(info_approximation))
361 WRITE (unit_nr_file,
'(A20,1X,10(2X,A20,1X))')
"# Frequency (eV)", αω
"Im{_{avg}()}", αω
"Im{_xx()}", &
362 αω
"Im{_xy()}", αω
"Im{_xz()}", αω
"Im{_yx()}", αω
"Im{_yy()}", αω
"Im{_yz()}", αω
"Im{_zx()}", &
363 αω
"Im{_zy()}", αω
"Im{_zz()}"
365 WRITE (unit_nr_file,
'(11(F20.8,1X))') abs_spectrum(i, 1)*
evolt, abs_spectrum(i, 2:11)
369 DEALLOCATE (abs_spectrum)
372 IF (unit_nr > 0)
THEN
373 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
374 WRITE (unit_nr,
'(T2,A4,T7,A,A)') &
375 'BSE|',
"Printed optical absorption spectrum to local files, e.g. "
376 WRITE (unit_nr,
'(T2,A4,T7,A)') &
377 'BSE|', file_name_spectrum
378 WRITE (unit_nr,
'(T2,A4,T7,A,A)') &
379 'BSE|',
"as well as photoabsorption cross section to, e.g. "
380 WRITE (unit_nr,
'(T2,A4,T7,A)') &
381 'BSE|', file_name_crosssection
382 WRITE (unit_nr,
'(T2,A4,T7,A52)') &
383 'BSE|',
"using the Eq. (7.51) from C. Ullrichs Book on TDDFT:"
384 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
385 WRITE (unit_nr,
'(T2,A4,T10,A75)') &
386 'BSE|', αωωη²Ω²
"Im{_{avg}()} = -Im{\sum_{n=1}^{N_exc} \frac{f_n}{(+i) - (^n)}}"
387 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
388 WRITE (unit_nr,
'(T2,A4,T7,A)') &
389 'BSE|',
"or for the full polarizability tensor:"
390 WRITE (unit_nr,
'(T2,A4,T10,A)') &
391 'BSE|', αµµωΩµµωη²Ω²
"_{ '}() = -\sum_n [2 ^n d_^n d_'^n] / [(+i)- (^n)]"
392 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
393 WRITE (unit_nr,
'(T2,A4,T7,A)') &
394 'BSE|',
"as well as Eq. (7.48):"
395 WRITE (unit_nr,
'(T2,A4,T10,A)') &
396 'BSE|', σµµωπωαµµω
"_{ '}() = 4 Im{_{ '}()} / c"
397 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
398 WRITE (unit_nr,
'(T2,A4,T7,A)') &
399 'BSE|', µ
"with transition moments d_^n, oscillator strengths f_n,"
400 WRITE (unit_nr,
'(T2,A4,T7,A)') &
401 'BSE|', Ω
"excitation energies ^n and the speed of light c."
402 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
403 WRITE (unit_nr,
'(T2,A4,T7,A)') &
404 'BSE|',
"Please note that we adopt an additional minus sign for both quantities,"
405 WRITE (unit_nr,
'(T2,A4,T7,A)') &
406 'BSE|',
"due to the convention for charge vs particle density as done in MolGW:"
407 WRITE (unit_nr,
'(T2,A4,T7,A)') &
408 'BSE|',
"https://doi.org/10.1016/j.cpc.2016.06.019."
409 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
412 CALL timestop(handle)
430 mo_coeff, homo, virtual, &
431 info_approximation, &
433 qs_env, unit_nr, mp2_env)
436 TYPE(
cp_fm_type),
INTENT(IN),
OPTIONAL :: fm_y
437 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(IN) :: mo_coeff
438 INTEGER,
INTENT(IN) :: homo, virtual
439 CHARACTER(LEN=10) :: info_approximation
440 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: oscill_str
442 INTEGER,
INTENT(IN) :: unit_nr
443 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
445 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_NTOs'
447 CHARACTER(LEN=20),
DIMENSION(2) :: nto_name
448 INTEGER :: handle, homo_irred, i, i_nto, info_svd, &
449 j, n_exc, n_nto, nao_full, nao_trunc
450 INTEGER,
DIMENSION(:),
POINTER :: stride
451 LOGICAL :: append_cube, cube_file
452 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eigval_svd_squ
453 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigval_svd
455 fm_struct_nto_holes, &
456 fm_struct_nto_particles, &
458 TYPE(
cp_fm_type) :: fm_eigvl, fm_eigvr_t, fm_m, fm_mo_coeff, fm_nto_coeff_holes, &
459 fm_nto_coeff_particles, fm_nto_set, fm_x_ia, fm_y_ai
460 TYPE(
mo_set_type),
ALLOCATABLE,
DIMENSION(:) :: nto_set
463 CALL timeset(routinen, handle)
468 nao_full = qs_env%mos(1)%nao
469 nao_trunc = homo + virtual
471 homo_irred = qs_env%mos(1)%homo
480 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
481 nao_trunc, nao_trunc)
483 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
486 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
489 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
492 CALL cp_fm_create(fm_mo_coeff, matrix_struct=fm_struct_mo_coeff, &
498 nao_full, nao_trunc, &
499 1, homo_irred - homo + 1, &
501 mo_coeff(1)%matrix_struct%context)
504 IF (unit_nr > 0)
THEN
505 WRITE (unit_nr,
'(T2,A4,T7,A)')
'BSE|', &
506 φχ
'The Natural Transition Orbital (NTO) pairs _I(r_e) and _I(r_h) for a fixed'
507 WRITE (unit_nr,
'(T2,A4,T7,A)')
'BSE|', &
508 'excitation index n are obtained by singular value decomposition of T'
509 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
510 WRITE (unit_nr,
'(T2,A4,T15,A)')
'BSE|', &
512 IF (
PRESENT(fm_y))
THEN
513 WRITE (unit_nr,
'(T2,A4,T15,A)')
'BSE|', &
516 WRITE (unit_nr,
'(T2,A4,T15,A)')
'BSE|', &
519 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
520 WRITE (unit_nr,
'(T2,A4,T15,A)')
'BSE|', &
522 WRITE (unit_nr,
'(T2,A4,T15,A)')
'BSE|', &
523 φψ
'_I(r_e) = \sum_p V_pI _p(r_e)'
524 WRITE (unit_nr,
'(T2,A4,T15,A)')
'BSE|', &
525 χψ
'_I(r_h) = \sum_p U_pI _p(r_e)'
526 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
527 WRITE (unit_nr,
'(T2,A4,T7,A)')
'BSE|', &
528 'where we have introduced'
529 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
530 WRITE (unit_nr,
'(T2,A4,T7,A,T20,A)') &
531 'BSE|', ψ
"_p:",
"occupied and virtual molecular orbitals,"
532 WRITE (unit_nr,
'(T2,A4,T7,A,T20,A)') &
533 'BSE|', φ
"_I(r_e):",
"NTO state for the electron,"
534 WRITE (unit_nr,
'(T2,A4,T7,A,T20,A)') &
535 'BSE|', χ
"_I(r_h):",
"NTO state for the hole,"
536 WRITE (unit_nr,
'(T2,A4,T7,A,T20,A)') &
537 'BSE|', Λ
":", λ
"diagonal matrix of NTO weights _I,"
538 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
539 WRITE (unit_nr,
'(T2,A4,T7,A)')
'BSE|', &
540 "The NTOs are calculated with the following settings:"
541 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
542 WRITE (unit_nr,
'(T2,A4,T7,A,T71,I10)')
'BSE|',
'Number of excitations, for which NTOs are computed', &
543 mp2_env%bse%num_print_exc_ntos
544 IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp)
THEN
545 WRITE (unit_nr,
'(T2,A4,T7,A,T71,F10.3)')
'BSE|',
'Threshold for oscillator strength f^n', &
546 mp2_env%bse%eps_nto_osc_str
548 WRITE (unit_nr,
'(T2,A4,T7,A,T71,A10)')
'BSE|',
'Threshold for oscillator strength f^n', &
551 WRITE (unit_nr,
'(T2,A4,T7,A,T72,F10.3)')
'BSE|', λ
'Threshold for NTO weights (_I)^2', &
552 mp2_env%bse%eps_nto_eigval
556 IF (unit_nr > 0)
THEN
557 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
558 IF (.NOT.
PRESENT(fm_y))
THEN
559 WRITE (unit_nr,
'(T2,A4,T7,A)')
'BSE|', &
560 'NTOs from solving the BSE within the TDA:'
562 WRITE (unit_nr,
'(T2,A4,T7,A)')
'BSE|', &
563 'NTOs from solving the BSE without the TDA:'
565 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
566 WRITE (unit_nr,
'(T2,A4,T8,A12,T22,A8,T33,A14,T62,A)')
'BSE|', &
567 'Excitation n',
"TDA/ABBA",
"Index of NTO I", λ
'NTO weights (_I)^2'
570 DO j = 1, mp2_env%bse%num_print_exc_ntos
571 n_exc = mp2_env%bse%bse_nto_state_list_final(j)
573 IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp)
THEN
575 IF (oscill_str(n_exc) < mp2_env%bse%eps_nto_osc_str)
THEN
577 IF (unit_nr > 0)
THEN
578 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
579 WRITE (unit_nr,
'(T2,A4,T8,I12,T24,A6,T42,A39)')
'BSE|', &
580 n_exc, info_approximation,
"Skipped (Oscillator strength too small)"
587 name=
"single_part_trans_dm")
590 CALL cp_fm_create(fm_nto_coeff_holes, matrix_struct=fm_struct_nto_holes, &
591 name=
"nto_coeffs_holes")
594 CALL cp_fm_create(fm_nto_coeff_particles, matrix_struct=fm_struct_nto_particles, &
595 name=
"nto_coeffs_particles")
600 .false., unit_nr, mp2_env)
609 IF (
PRESENT(fm_y))
THEN
612 .true., unit_nr, mp2_env)
628 matrix_struct=fm_m%matrix_struct, &
629 name=
"LEFT_SINGULAR_MATRIX")
632 matrix_struct=fm_m%matrix_struct, &
633 name=
"RIGHT_SINGULAR_MATRIX")
636 ALLOCATE (eigval_svd(nao_trunc))
637 eigval_svd(:) = 0.0_dp
639 CALL cp_fm_svd(fm_m, fm_eigvl, fm_eigvr_t, eigval_svd, info_svd)
640 IF (info_svd /= 0)
THEN
641 IF (unit_nr > 0)
THEN
642 CALL cp_warn(__location__, &
643 "SVD for computation of NTOs not successful. "// &
644 "Skipping print of NTOs.")
645 IF (info_svd > 0)
THEN
646 CALL cp_warn(__location__, &
647 "PDGESVD detected heterogeneity. "// &
648 "Decreasing number of MPI ranks might solve this issue.")
657 ALLOCATE (eigval_svd_squ(nao_trunc))
658 eigval_svd_squ(:) = eigval_svd(:)**2
660 IF (.NOT.
PRESENT(fm_y))
THEN
661 IF (abs(sum(eigval_svd_squ) - 1) >= 1e-5)
THEN
662 cpwarn(
"Sum of NTO coefficients deviates from 1!")
668 CALL parallel_gemm(
"N",
"N", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvl, 0.0_dp, &
672 CALL parallel_gemm(
"N",
"T", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvr_t, 0.0_dp, &
673 fm_nto_coeff_particles)
681 nto_name(1) =
'Hole_coord'
682 nto_name(2) =
'Particle_coord'
683 ALLOCATE (nto_set(2))
686 DO i_nto = 1, nao_trunc
687 IF (eigval_svd_squ(i_nto) > mp2_env%bse%eps_nto_eigval)
THEN
695 IF (unit_nr > 0)
THEN
696 WRITE (unit_nr,
'(T2,A4)')
'BSE|'
698 WRITE (unit_nr,
'(T2,A4,T8,I12,T24,A6,T41,I6,T71,F10.5)')
'BSE|', &
699 n_exc, info_approximation, i_nto, eigval_svd_squ(i_nto)
707 CALL allocate_mo_set(nto_set(i), nao_trunc, n_nto, 0, 0.0_dp, 2.0_dp, 0.0_dp)
708 CALL init_mo_set(nto_set(i), fm_ref=fm_nto_set, name=nto_name(i))
714 CALL cp_fm_to_fm(fm_nto_coeff_holes, nto_set(1)%mo_coeff, ncol=n_nto)
715 CALL cp_fm_to_fm(fm_nto_coeff_particles, nto_set(2)%mo_coeff, ncol=n_nto)
724 stride, append_cube, nto_section)
729 DEALLOCATE (eigval_svd)
730 DEALLOCATE (eigval_svd_squ)
744 CALL timestop(handle)
761 fm_multipole_ij_trunc, fm_multipole_ab_trunc, &
762 fm_multipole_ai_trunc, &
763 i_exc, homo, virtual, &
767 DIMENSION(:) :: exc_descr
769 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
770 INTENT(IN) :: fm_multipole_ij_trunc, &
771 fm_multipole_ab_trunc, &
772 fm_multipole_ai_trunc
773 INTEGER,
INTENT(IN) :: i_exc, homo, virtual
774 TYPE(
cp_fm_type),
INTENT(IN),
OPTIONAL :: fm_y_ia
776 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_exciton_descriptors'
778 INTEGER :: handle, i_dir, j_dir
779 INTEGER,
DIMENSION(3) :: mask_quadrupole
781 REAL(kind=
dp) :: norm_x, norm_xpy, norm_y
782 REAL(kind=
dp),
DIMENSION(3) :: r_e_sq_x, r_e_sq_y, r_e_x, r_e_y, &
783 r_h_sq_x, r_h_sq_y, r_h_x, r_h_y
784 REAL(kind=
dp),
DIMENSION(3, 3) :: r_e_h_xx, r_e_h_xy, r_e_h_yx, r_e_h_yy
786 TYPE(
cp_fm_type) :: fm_work_ba, fm_work_ia, fm_work_ia_2
788 CALL timeset(routinen, handle)
789 IF (
PRESENT(fm_y_ia))
THEN
797 mask_quadrupole = [4, 7, 9]
800 context=fm_x_ia%matrix_struct%context, nrow_global=homo, ncol_global=virtual)
802 context=fm_x_ia%matrix_struct%context, nrow_global=virtual, ncol_global=virtual)
812 r_e_h_xx(:, :) = 0.0_dp
813 r_e_h_xy(:, :) = 0.0_dp
814 r_e_h_yx(:, :) = 0.0_dp
815 r_e_h_yy(:, :) = 0.0_dp
822 exc_descr(i_exc)%r_e(:) = 0.0_dp
823 exc_descr(i_exc)%r_h(:) = 0.0_dp
824 exc_descr(i_exc)%r_e_sq(:) = 0.0_dp
825 exc_descr(i_exc)%r_h_sq(:) = 0.0_dp
826 exc_descr(i_exc)%r_e_h(:, :) = 0.0_dp
828 exc_descr(i_exc)%flag_TDA = flag_tda
829 exc_descr(i_exc)%norm_XpY = 0.0_dp
835 IF (.NOT. flag_tda)
THEN
837 norm_xpy = norm_xpy + norm_y
840 exc_descr(i_exc)%norm_XpY = norm_xpy
846 r_h_x(i_dir) = r_h_x(i_dir)/norm_xpy
847 IF (.NOT. flag_tda)
THEN
850 r_h_y(i_dir) = r_h_y(i_dir)/norm_xpy
853 exc_descr(i_exc)%r_h(:) = r_h_x(:) + r_h_y(:)
859 r_e_x(i_dir) = r_e_x(i_dir)/norm_xpy
860 IF (.NOT. flag_tda)
THEN
863 r_e_y(i_dir) = r_e_y(i_dir)/norm_xpy
866 exc_descr(i_exc)%r_e(:) = r_e_x(:) + r_e_y(:)
872 fm_x_ia, r_h_sq_x(i_dir))
873 r_h_sq_x(i_dir) = r_h_sq_x(i_dir)/norm_xpy
874 IF (.NOT. flag_tda)
THEN
877 fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_h_sq_y(i_dir))
878 r_h_sq_y(i_dir) = r_h_sq_y(i_dir)/norm_xpy
881 exc_descr(i_exc)%r_h_sq(:) = r_h_sq_x(:) + r_h_sq_y(:)
887 fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_e_sq_x(i_dir))
888 r_e_sq_x(i_dir) = r_e_sq_x(i_dir)/norm_xpy
889 IF (.NOT. flag_tda)
THEN
892 fm_y_ia, r_e_sq_y(i_dir))
893 r_e_sq_y(i_dir) = r_e_sq_y(i_dir)/norm_xpy
896 exc_descr(i_exc)%r_e_sq(:) = r_e_sq_x(:) + r_e_sq_y(:)
911 CALL parallel_gemm(
"N",
"N", homo, virtual, virtual, 1.0_dp, &
912 fm_x_ia, fm_multipole_ab_trunc(i_dir), 0.0_dp, fm_work_ia)
915 fm_multipole_ij_trunc(j_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
917 CALL cp_fm_trace(fm_x_ia, fm_work_ia_2, r_e_h_xx(i_dir, j_dir))
918 r_e_h_xx(i_dir, j_dir) = r_e_h_xx(i_dir, j_dir)/norm_xpy
919 IF (.NOT. flag_tda)
THEN
924 CALL parallel_gemm(
"N",
"N", homo, virtual, virtual, 1.0_dp, &
925 fm_y_ia, fm_multipole_ab_trunc(j_dir), 0.0_dp, fm_work_ia)
928 fm_multipole_ij_trunc(i_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
930 CALL cp_fm_trace(fm_y_ia, fm_work_ia_2, r_e_h_yy(i_dir, j_dir))
931 r_e_h_yy(i_dir, j_dir) = r_e_h_yy(i_dir, j_dir)/norm_xpy
945 CALL parallel_gemm(
"T",
"T", virtual, virtual, homo, 1.0_dp, &
946 fm_y_ia, fm_multipole_ai_trunc(j_dir), 0.0_dp, fm_work_ba)
948 CALL parallel_gemm(
"T",
"N", homo, virtual, virtual, 1.0_dp, &
949 fm_multipole_ai_trunc(i_dir), fm_work_ba, 0.0_dp, fm_work_ia)
951 CALL cp_fm_trace(fm_work_ia, fm_x_ia, r_e_h_xy(i_dir, j_dir))
952 r_e_h_xy(i_dir, j_dir) = r_e_h_xy(i_dir, j_dir)/norm_xpy
966 CALL parallel_gemm(
"T",
"T", virtual, virtual, homo, 1.0_dp, &
967 fm_x_ia, fm_multipole_ai_trunc(i_dir), 0.0_dp, fm_work_ba)
969 CALL parallel_gemm(
"T",
"N", homo, virtual, virtual, 1.0_dp, &
970 fm_multipole_ai_trunc(j_dir), fm_work_ba, 0.0_dp, fm_work_ia)
972 CALL cp_fm_trace(fm_work_ia, fm_y_ia, r_e_h_yx(i_dir, j_dir))
973 r_e_h_yx(i_dir, j_dir) = r_e_h_yx(i_dir, j_dir)/norm_xpy
977 exc_descr(i_exc)%r_e_h(:, :) = r_e_h_xx(:, :) + r_e_h_xy(:, :) + r_e_h_yx(:, :) + r_e_h_yy(:, :)
987 exc_descr(i_exc)%diff_r_abs = sqrt(sum((exc_descr(i_exc)%r_h(:) - exc_descr(i_exc)%r_e(:))**2))
990 exc_descr(i_exc)%sigma_e = sqrt(sum(exc_descr(i_exc)%r_e_sq(:)) - sum(exc_descr(i_exc)%r_e(:)**2))
993 exc_descr(i_exc)%sigma_h = sqrt(sum(exc_descr(i_exc)%r_h_sq(:)) - sum(exc_descr(i_exc)%r_h(:)**2))
997 exc_descr(i_exc)%d_eh_dir(i_dir) = abs(exc_descr(i_exc)%r_h(i_dir) - exc_descr(i_exc)%r_e(i_dir))
998 exc_descr(i_exc)%sigma_e_dir(i_dir) = sqrt(exc_descr(i_exc)%r_e_sq(i_dir) - exc_descr(i_exc)%r_e(i_dir)**2)
999 exc_descr(i_exc)%sigma_h_dir(i_dir) = sqrt(exc_descr(i_exc)%r_h_sq(i_dir) - exc_descr(i_exc)%r_h(i_dir)**2)
1004 exc_descr(i_exc)%cov_e_h_sum = 0.0_dp
1005 exc_descr(i_exc)%cov_e_h(:, :) = 0.0_dp
1006 exc_descr(i_exc)%corr_e_h_matrix(:, :) = 0.0_dp
1009 exc_descr(i_exc)%cov_e_h(i_dir, j_dir) = exc_descr(i_exc)%r_e_h(i_dir, j_dir) &
1010 - exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(j_dir)
1011 exc_descr(i_exc)%corr_e_h_matrix(i_dir, j_dir) = &
1012 exc_descr(i_exc)%cov_e_h(i_dir, j_dir)/ &
1013 (exc_descr(i_exc)%sigma_e_dir(i_dir)*exc_descr(i_exc)%sigma_h_dir(j_dir))
1015 exc_descr(i_exc)%cov_e_h_sum = exc_descr(i_exc)%cov_e_h_sum + &
1016 exc_descr(i_exc)%r_e_h(i_dir, i_dir) - &
1017 exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(i_dir)
1021 exc_descr(i_exc)%corr_e_h = exc_descr(i_exc)%cov_e_h_sum/(exc_descr(i_exc)%sigma_e*exc_descr(i_exc)%sigma_h)
1024 exc_descr(i_exc)%diff_r_sqr = sqrt(exc_descr(i_exc)%diff_r_abs**2 + &
1025 exc_descr(i_exc)%sigma_e**2 + exc_descr(i_exc)%sigma_h**2 &
1026 - 2*exc_descr(i_exc)%cov_e_h_sum)
1029 exc_descr(i_exc)%d_exc_dir(i_dir) = sqrt(exc_descr(i_exc)%d_eh_dir(i_dir)**2 + &
1030 exc_descr(i_exc)%sigma_e_dir(i_dir)**2 + &
1031 exc_descr(i_exc)%sigma_h_dir(i_dir)**2 - &
1032 2*exc_descr(i_exc)%cov_e_h(i_dir, i_dir))
1036 exc_descr(i_exc)%r_e_shift(:) = exc_descr(i_exc)%r_e(:)
1037 exc_descr(i_exc)%r_h_shift(:) = exc_descr(i_exc)%r_h(:)
1042 CALL timestop(handle)
Routines for computing excitonic properties, e.g. exciton diameter, from the BSE.
subroutine, public compute_and_print_absorption_spectrum(oscill_str, polarizability_residues, exc_ens, info_approximation, unit_nr, mp2_env)
Computes and returns absorption spectrum for the frequency range and broadening provided by the user....
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 trace_exciton_descr(fm_a, fm_b, fm_c, alpha)
Computes trace of form Tr{A^T B C} for exciton descriptors.
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, unit_nr, mp2_env)
...
subroutine, public print_bse_nto_cubes(qs_env, mos, istate, info_approximation, stride, append_cube, print_section)
Borrowed from the tddfpt module with slight adaptions.
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)
Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,...
Utility routines to open and close files. Tracking of preconnections.
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
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_svd(matrix_a, matrix_eigvl, matrix_eigvr_t, eigval, info)
decomposes a quadratic matrix into its singular value decomposition
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_create(matrix, matrix_struct, name, use_sp, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_to_fm_submat_general(source, destination, nrows, ncols, s_firstrow, s_firstcol, d_firstrow, d_firstcol, global_context)
General copy of a submatrix of fm matrix to a submatrix of another fm matrix. The two matrices can ha...
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_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,...
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Types needed for MP2 calculations.
basic linear algebra operations for full matrixes
Definition of physical constants:
real(kind=dp), parameter, public c_light_au
real(kind=dp), parameter, public evolt
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, 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, 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, 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.
Definition and initialisation of the mo data type.
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...