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_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,...
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
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...