64#include "./base/base_uses.f90"
70 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_tddfpt2_eigensolver'
72 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
74 INTEGER,
PARAMETER,
PRIVATE :: nderivs = 3
75 INTEGER,
PARAMETER,
PRIVATE :: maxspins = 2
106 tddfpt_control, S_C0)
107 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: evects
108 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(in) :: s_c0_c0t
112 REAL(kind=
dp),
DIMENSION(:),
INTENT(in) :: evals
114 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(in) :: s_c0
116 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_orthogonalize_psi1_psi0'
118 INTEGER :: handle, ispin, ivect, nactive, nao, &
122 TYPE(
cp_fm_type),
DIMENSION(:),
POINTER :: mos
125 CALL timeset(routinen, handle)
127 nspins =
SIZE(evects, 1)
128 nvects =
SIZE(evects, 2)
131 IF (.NOT. tddfpt_control%do_smearing)
THEN
138 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
139 nrow_global=nao, ncol_global=nactive)
143 CALL parallel_gemm(
'T',
'N', nao, nactive, nao, 1.0_dp, s_c0_c0t(spin), &
144 evects(ispin, ivect), 0.0_dp, evortho)
153 ALLOCATE (mos(nspins))
155 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
156 nrow_global=nao, ncol_global=nactive)
171 CALL timestop(handle)
189 FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm, spinflip) &
191 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: evects
192 TYPE(
cp_fm_type),
DIMENSION(:),
INTENT(in) :: s_c0
193 REAL(kind=
dp),
INTENT(in) :: max_norm
194 INTEGER,
INTENT(in) :: spinflip
195 LOGICAL :: is_nonortho
197 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_is_nonorthogonal_psi1_psi0'
199 INTEGER :: handle, ispin, ivect, nactive, nao, &
200 nocc, nspins, nvects, spin
201 REAL(kind=
dp) :: maxabs_val
205 CALL timeset(routinen, handle)
207 nspins =
SIZE(evects, 1)
208 nvects =
SIZE(evects, 2)
210 is_nonortho = .false.
212 loop:
DO ispin = 1, nspins
221 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
222 nrow_global=nao, ncol_global=nactive)
224 ncol_global=nactive, template_fmstruct=matrix_struct)
229 CALL parallel_gemm(
'T',
'N', nocc, nactive, nao, 1.0_dp, s_c0(spin), &
230 evects(ispin, ivect), 0.0_dp, aortho)
232 is_nonortho = maxabs_val > max_norm
233 IF (is_nonortho)
THEN
241 CALL timestop(handle)
243 END FUNCTION tddfpt_is_nonorthogonal_psi1_psi0
280 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: evects
281 INTEGER,
INTENT(in) :: nvects_new
282 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(INOUT) :: s_evects
285 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_orthonormalize_psi1_psi1'
287 INTEGER :: handle, ispin, ivect, jvect, nspins, &
288 nvects_old, nvects_total
289 INTEGER,
DIMENSION(maxspins) :: nactive
290 REAL(kind=
dp) :: norm
291 REAL(kind=
dp),
DIMENSION(maxspins) :: weights
293 CALL timeset(routinen, handle)
295 nspins =
SIZE(evects, 1)
296 nvects_total =
SIZE(evects, 2)
297 nvects_old = nvects_total - nvects_new
299 IF (debug_this_module)
THEN
300 cpassert(
SIZE(s_evects, 1) == nspins)
301 cpassert(
SIZE(s_evects, 2) == nvects_total)
302 cpassert(nvects_old >= 0)
306 CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
309 DO jvect = nvects_old + 1, nvects_total
311 DO ivect = 1, jvect - 1
312 CALL cp_fm_trace(evects(:, jvect), s_evects(:, ivect), weights(1:nspins), accurate=.false.)
313 norm = sum(weights(1:nspins))
323 ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
326 CALL cp_fm_trace(evects(:, jvect), s_evects(:, jvect), weights(1:nspins), accurate=.false.)
328 norm = sum(weights(1:nspins))
329 norm = 1.0_dp/sqrt(norm)
337 CALL timestop(handle)
359 SUBROUTINE tddfpt_compute_aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, &
360 matrix_ks, qs_env, kernel_env, &
361 sub_env, work_matrices, matrix_s)
362 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(INOUT) :: aop_evects
363 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: evects, s_evects
367 TYPE(
dbcsr_p_type),
DIMENSION(:),
INTENT(in) :: matrix_ks
374 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_Aop_evects'
376 INTEGER :: handle, ispin, ivect, nspins, nvects
377 INTEGER,
DIMENSION(maxspins) :: nmo_occ
378 LOGICAL :: do_admm, do_bse, do_bse_gw_only, &
379 do_bse_w_only, do_hfx, &
380 do_lri_response, is_rks_triplets, &
382 REAL(kind=
dp) :: rcut, scale
387 CALL timeset(routinen, handle)
389 nspins =
SIZE(gs_mos, 1)
390 nvects =
SIZE(evects, 2)
391 do_hfx = tddfpt_control%do_hfx
392 do_admm = tddfpt_control%do_admm
394 kernel_env_admm_aux => kernel_env%admm_kernel
396 NULLIFY (kernel_env_admm_aux)
398 is_rks_triplets = tddfpt_control%rks_triplets
399 do_lri_response = tddfpt_control%do_lrigpw
400 do_bse = tddfpt_control%do_bse
401 do_bse_w_only = tddfpt_control%do_bse_w_only
402 do_bse_gw_only = tddfpt_control%do_bse_gw_only
403 IF (do_bse) do_hfx = .false.
404 IF (do_bse_w_only) do_hfx = .false.
405 IF (do_bse_gw_only) do_hfx = .false.
407 IF (debug_this_module)
THEN
409 cpassert(
SIZE(aop_evects, 1) ==
SIZE(evects, 1))
410 cpassert(
SIZE(s_evects, 1) ==
SIZE(evects, 1))
411 cpassert(
SIZE(aop_evects, 2) == nvects)
412 cpassert(
SIZE(s_evects, 2) == nvects)
413 cpassert(
SIZE(gs_mos) == nspins)
417 nmo_occ(ispin) =
SIZE(gs_mos(ispin)%evals_occ)
422 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
424 DO ispin = 1,
SIZE(evects, 1)
425 associate(evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect))
426 IF (
ASSOCIATED(evect%matrix_struct))
THEN
427 IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
432 ELSE IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
444 CALL fhxc_kernel(aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, &
445 kernel_env%full_kernel, kernel_env_admm_aux, sub_env, work_matrices, &
446 tddfpt_control%admm_symm, tddfpt_control%admm_xc_correction, &
447 do_lri_response, tddfpt_mgrid=tddfpt_control%mgrid_is_explicit)
450 CALL stda_kernel(aop_evects, evects, is_rks_triplets, qs_env, tddfpt_control%stda_control, &
451 kernel_env%stda_kernel, sub_env, work_matrices)
455 DO ispin = 1,
SIZE(evects, 1)
460 cpabort(
"Kernel type not implemented")
463 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
465 DO ispin = 1,
SIZE(evects, 1)
466 associate(aop_evect => aop_evects(ispin, ivect), &
467 work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
468 IF (
ASSOCIATED(aop_evect%matrix_struct))
THEN
469 IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
474 ELSE IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
485 IF ((.NOT. do_bse) .AND. (.NOT. do_bse_gw_only))
THEN
486 CALL tddfpt_apply_energy_diff(aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
487 gs_mos=gs_mos, matrix_ks=matrix_ks, tddfpt_control=tddfpt_control)
489 CALL zeroth_order_gw(qs_env=qs_env, aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
490 gs_mos=gs_mos, matrix_s=matrix_s, matrix_ks=matrix_ks)
494 IF (tddfpt_control%do_smearing)
THEN
497 CALL add_smearing_aterm(qs_env, aop_evects(ispin, ivect), evects(ispin, ivect), &
498 s_evects(ispin, ivect), gs_mos(ispin)%mos_occ, &
499 tddfpt_control%smeared_occup(ispin)%fermia, matrix_s)
505 IF (tddfpt_control%kernel == tddfpt_kernel_full)
THEN
507 CALL tddfpt_apply_hfx(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
508 qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
509 work_hmat_symm=work_matrices%hfx_hmat_symm, &
510 work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
511 work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
512 work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
513 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda)
THEN
516 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none)
THEN
520 cpabort(
"Kernel type not implemented")
524 IF (tddfpt_control%kernel == tddfpt_kernel_full)
THEN
525 IF (tddfpt_control%do_hfxsr)
THEN
526 re_int = tddfpt_control%hfxsr_re_int
528 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
529 kernel_env%full_kernel%admm_env, &
530 kernel_env%full_kernel%hfxsr_section, &
531 kernel_env%full_kernel%x_data, 1, re_int, &
532 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
533 work_hmat=work_matrices%hfxsr_hmat_symm, &
534 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
536 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
537 kernel_env%full_kernel%admm_env, &
538 kernel_env%full_kernel%hfxsr_section, &
539 kernel_env%full_kernel%x_data, -1, .false., &
540 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
541 work_hmat=work_matrices%hfxsr_hmat_asymm, &
542 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
543 tddfpt_control%hfxsr_re_int = .false.
545 IF (tddfpt_control%do_hfxlr)
THEN
546 rcut = tddfpt_control%hfxlr_rcut
547 scale = tddfpt_control%hfxlr_scale
549 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
550 IF (
ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct))
THEN
551 CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
552 work_matrices%evects_sub(:, ivect), &
553 work_matrices%Aop_evects_sub(:, ivect))
559 CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
560 evects(:, ivect), aop_evects(:, ivect))
565 IF ((do_bse) .OR. (do_bse_w_only))
THEN
568 CALL tddfpt_apply_bse(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env, &
574 CALL timestop(handle)
576 END SUBROUTINE tddfpt_compute_aop_evects
593 SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
595 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: ritz_vects, aop_ritz
596 REAL(kind=dp),
DIMENSION(:),
INTENT(out) :: evals
597 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: krylov_vects, aop_krylov
598 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: atilde
599 INTEGER,
INTENT(IN) :: nkvo, nkvn
601 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_ritz_vects'
603 INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, &
606 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: at12, at21, at22, evects_atilde
607 TYPE(cp_blacs_env_type),
POINTER :: blacs_env_global
608 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
609 TYPE(cp_fm_type) :: atilde_fm, evects_atilde_fm
611 CALL timeset(routinen, handle)
613 nspins =
SIZE(krylov_vects, 1)
614 nkvs =
SIZE(krylov_vects, 2)
615 nrvs =
SIZE(ritz_vects, 2)
616 cpassert(nkvs == nkvo + nkvn)
618 CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
620 CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
621 CALL cp_fm_create(atilde_fm, fm_struct, set_zero=.true.)
622 CALL cp_fm_create(evects_atilde_fm, fm_struct, set_zero=.true.)
623 CALL cp_fm_struct_release(fm_struct)
626 CALL reallocate(atilde, 1, nkvs, 1, nkvs)
631 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
632 atilde(1:nkvs, 1:nkvs), accurate=.false.)
634 ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
635 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
636 at12, accurate=.false.)
637 atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
638 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
639 at21, accurate=.false.)
640 atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
641 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
642 at22, accurate=.false.)
643 atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
644 DEALLOCATE (at12, at21, at22)
646 atilde = 0.5_dp*(atilde + transpose(atilde))
647 CALL cp_fm_set_submatrix(atilde_fm, atilde)
650 CALL choose_eigv_solver(atilde_fm, evects_atilde_fm, evals(1:nkvs))
652 ALLOCATE (evects_atilde(nkvs, nrvs))
653 CALL cp_fm_get_submatrix(evects_atilde_fm, evects_atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
654 CALL cp_fm_release(evects_atilde_fm)
655 CALL cp_fm_release(atilde_fm)
662 CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
663 CALL cp_fm_set_all(aop_ritz(ispin, irv), 0.0_dp)
667 act = evects_atilde(ikv, irv)
669 CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
670 act, krylov_vects(ispin, ikv))
671 CALL cp_fm_scale_and_add(1.0_dp, aop_ritz(ispin, irv), &
672 act, aop_krylov(ispin, ikv))
678 DEALLOCATE (evects_atilde)
680 CALL timestop(handle)
682 END SUBROUTINE tddfpt_compute_ritz_vects
697 SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
698 matrix_s, tddfpt_control)
699 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: residual_vects
700 REAL(kind=dp),
DIMENSION(:),
INTENT(in) :: evals
701 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: ritz_vects, aop_ritz
702 TYPE(tddfpt_ground_state_mos),
DIMENSION(:), &
704 TYPE(dbcsr_type),
POINTER :: matrix_s
705 TYPE(tddfpt2_control_type),
POINTER :: tddfpt_control
707 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_residual_vects'
708 REAL(kind=dp),
PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*epsilon(1.0_dp)
710 INTEGER :: handle, ica, icb, icol_local, &
711 irow_local, irv, ispin, nao, &
712 ncols_local, nrows_local, nrvs, &
713 nspins, spin2, spinflip
714 INTEGER,
DIMENSION(:),
POINTER :: col_indices_local, row_indices_local
715 INTEGER,
DIMENSION(maxspins) :: nactive, nmo_virt
716 REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
717 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :), &
718 POINTER :: weights_ldata
719 TYPE(cp_fm_struct_type),
POINTER :: ao_mo_struct, virt_mo_struct
720 TYPE(cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: awork, vomat
722 CALL timeset(routinen, handle)
724 nspins =
SIZE(residual_vects, 1)
725 nrvs =
SIZE(residual_vects, 2)
726 spinflip = tddfpt_control%spinflip
729 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
730 ALLOCATE (awork(nspins), vomat(nspins))
732 IF (spinflip == no_sf_tddfpt)
THEN
737 nmo_virt(spin2) =
SIZE(gs_mos(spin2)%evals_virt)
739 CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
740 ncol_global=nactive(ispin))
741 CALL cp_fm_create(awork(ispin), ao_mo_struct)
743 CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(spin2), &
744 ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
745 CALL cp_fm_create(vomat(ispin), virt_mo_struct)
746 CALL cp_fm_struct_release(virt_mo_struct)
754 IF (spinflip == no_sf_tddfpt)
THEN
759 CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
760 ncol_local=ncols_local, row_indices=row_indices_local, &
761 col_indices=col_indices_local, local_data=weights_ldata)
764 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
765 ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
766 CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, aop_ritz(ispin, irv))
768 CALL parallel_gemm(
'T',
'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, gs_mos(spin2)%mos_virt, &
769 awork(ispin), 0.0_dp, vomat(ispin))
772 DO icol_local = 1, ncols_local
773 ica = col_indices_local(icol_local)
774 icb = gs_mos(ispin)%index_active(ica)
775 e_occ_plus_lambda = gs_mos(ispin)%evals_occ(icb) + lambda
777 DO irow_local = 1, nrows_local
778 eref = gs_mos(spin2)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
782 IF (abs(eref) < threshold) &
783 eref = eref + (1.0_dp - eref_scale)*lambda
785 weights_ldata(irow_local, icol_local) = weights_ldata(irow_local, icol_local)/eref
789 CALL parallel_gemm(
'N',
'N', nao, nactive(ispin), nmo_virt(spin2), 1.0_dp, gs_mos(spin2)%mos_virt, &
790 vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
794 CALL cp_fm_release(awork)
795 CALL cp_fm_release(vomat)
798 CALL timestop(handle)
800 END SUBROUTINE tddfpt_compute_residual_vects
826 matrix_ks, qs_env, kernel_env, &
827 sub_env, logger, iter_unit, energy_unit, &
828 tddfpt_print_section, work_matrices)
RESULT(conv)
829 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(inout) :: evects
830 REAL(kind=dp),
DIMENSION(:),
INTENT(inout) :: evals
831 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(inout) :: s_evects
832 TYPE(tddfpt_ground_state_mos),
DIMENSION(:), &
834 TYPE(tddfpt2_control_type),
POINTER :: tddfpt_control
835 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks
836 TYPE(qs_environment_type),
POINTER :: qs_env
837 TYPE(kernel_env_type),
INTENT(in) :: kernel_env
838 TYPE(tddfpt_subgroup_env_type),
INTENT(in) :: sub_env
839 TYPE(cp_logger_type),
POINTER :: logger
840 INTEGER,
INTENT(in) :: iter_unit, energy_unit
841 TYPE(section_vals_type),
POINTER :: tddfpt_print_section
842 TYPE(tddfpt_work_matrices),
INTENT(inout) :: work_matrices
843 REAL(kind=dp) :: conv
845 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_davidson_solver'
847 INTEGER :: handle, ispin, istate, iter, &
848 max_krylov_vects, nspins, nstates, &
849 nstates_conv, nvects_exist, nvects_new
850 INTEGER(kind=int_8) :: nstates_total
851 LOGICAL :: is_nonortho
852 REAL(kind=dp) :: t1, t2
853 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: evals_last
854 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: atilde
855 TYPE(cp_fm_type),
ALLOCATABLE,
DIMENSION(:, :) :: aop_krylov, aop_ritz, krylov_vects, &
857 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
859 CALL timeset(routinen, handle)
861 nspins =
SIZE(evects, 1)
862 nstates = tddfpt_control%nstates
863 nstates_total = tddfpt_total_number_of_states(tddfpt_control, gs_mos)
865 IF (debug_this_module)
THEN
866 cpassert(
SIZE(evects, 1) == nspins)
867 cpassert(
SIZE(evects, 2) == nstates)
868 cpassert(
SIZE(evals) == nstates)
871 CALL get_qs_env(qs_env, matrix_s=matrix_s)
874 max_krylov_vects = min(max(tddfpt_control%nkvs, nstates), int(nstates_total))
876 ALLOCATE (aop_ritz(nspins, nstates))
877 DO istate = 1, nstates
879 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_ritz(ispin, istate))
883 ALLOCATE (evals_last(max_krylov_vects))
884 ALLOCATE (aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), &
885 s_krylov(nspins, max_krylov_vects))
887 DO istate = 1, nstates
889 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
890 CALL cp_fm_to_fm(evects(ispin, istate), krylov_vects(ispin, istate))
892 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, s_krylov(ispin, istate))
893 CALL cp_fm_to_fm(s_evects(ispin, istate), s_krylov(ispin, istate))
895 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_krylov(ispin, istate))
904 ALLOCATE (atilde(1, 1))
908 CALL cp_iterate(logger%iter_info, iter_nr_out=iter)
911 CALL tddfpt_compute_aop_evects(aop_evects=aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
912 evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
913 s_evects=s_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
914 gs_mos=gs_mos, tddfpt_control=tddfpt_control, &
915 matrix_ks=matrix_ks, &
916 qs_env=qs_env, kernel_env=kernel_env, &
918 work_matrices=work_matrices, &
919 matrix_s=matrix_s(1)%matrix)
921 CALL tddfpt_compute_ritz_vects(ritz_vects=evects, aop_ritz=aop_ritz, &
922 evals=evals_last(1:nvects_exist + nvects_new), &
923 krylov_vects=krylov_vects(:, 1:nvects_exist + nvects_new), &
924 aop_krylov=aop_krylov(:, 1:nvects_exist + nvects_new), &
925 atilde=atilde, nkvo=nvects_exist, nkvn=nvects_new)
927 CALL tddfpt_write_restart(evects=evects, evals=evals_last(1:nstates), gs_mos=gs_mos, &
928 logger=logger, tddfpt_print_section=tddfpt_print_section)
930 conv = maxval(abs(evals_last(1:nstates) - evals(1:nstates)))
932 nvects_exist = nvects_exist + nvects_new
933 IF (nvects_exist + nvects_new > max_krylov_vects) &
934 nvects_new = max_krylov_vects - nvects_exist
935 IF (iter >= tddfpt_control%niters) nvects_new = 0
937 IF (conv > tddfpt_control%conv .AND. nvects_new > 0)
THEN
939 DO istate = 1, nvects_new
941 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
942 krylov_vects(ispin, nvects_exist + istate))
943 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
944 s_krylov(ispin, nvects_exist + istate))
945 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
946 aop_krylov(ispin, nvects_exist + istate))
950 CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
951 evals=evals_last(1:nvects_new), &
952 ritz_vects=evects(:, 1:nvects_new), aop_ritz=aop_ritz(:, 1:nvects_new), &
953 gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix, tddfpt_control=tddfpt_control)
956 work_matrices%S_C0_C0T, qs_env, &
957 gs_mos, evals(1:nstates), tddfpt_control, work_matrices%S_C0)
960 s_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
962 is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
963 work_matrices%S_C0, tddfpt_control%orthogonal_eps, &
964 tddfpt_control%spinflip)
968 is_nonortho = .false.
972 IF (energy_unit > 0)
THEN
973 WRITE (energy_unit,
'(/,4X,A,T14,A,T36,A)')
"State",
"Exc. energy (eV)",
"Convergence (eV)"
974 DO istate = 1, nstates
975 WRITE (energy_unit,
'(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
976 evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
978 WRITE (energy_unit, *)
979 CALL m_flush(energy_unit)
982 IF (iter_unit > 0)
THEN
984 DO istate = 1, nstates
985 IF (abs(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) &
986 nstates_conv = nstates_conv + 1
989 WRITE (iter_unit,
'(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
990 CALL m_flush(iter_unit)
994 evals(1:nstates) = evals_last(1:nstates)
997 IF (nvects_new == 0 .OR. is_nonortho)
THEN
1001 evals(1:nstates), tddfpt_control, work_matrices%S_C0)
1010 DO istate = nvects_exist + nvects_new, 1, -1
1011 DO ispin = nspins, 1, -1
1012 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_krylov(ispin, istate))
1013 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, s_krylov(ispin, istate))
1014 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
1017 DEALLOCATE (aop_krylov, krylov_vects, s_krylov)
1018 DEALLOCATE (evals_last)
1020 DO istate = nstates, 1, -1
1021 DO ispin = nspins, 1, -1
1022 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_ritz(ispin, istate))
1025 DEALLOCATE (aop_ritz)
1027 CALL timestop(handle)
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
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....
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
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...
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
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_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
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_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
subroutine, public cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
find the maximum absolute value of the matrix element maxval(abs(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 ...
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Utility routines for the memory handling.
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
Definition of physical constants:
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, 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, 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.
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
subroutine, public zeroth_order_gw(qs_env, aop_evects, evects, s_evects, gs_mos, matrix_s, matrix_ks)
...
subroutine, public tddfpt_apply_bse(aop_evects, evects, gs_mos, qs_env, s_evects)
...
subroutine, public tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, s_evects, matrix_s)
Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
subroutine, public tddfpt_orthogonalize_psi1_psi0(evects, s_c0_c0t, qs_env, gs_mos, evals, tddfpt_control, s_c0)
Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
real(kind=dp) function, public tddfpt_davidson_solver(evects, evals, s_evects, gs_mos, tddfpt_control, matrix_ks, qs_env, kernel_env, sub_env, logger, iter_unit, energy_unit, tddfpt_print_section, work_matrices)
Perform Davidson iterations.
subroutine, public fhxc_kernel(aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, kernel_env, kernel_env_admm_aux, sub_env, work_matrices, admm_symm, admm_xc_correction, do_lrigpw, tddfpt_mgrid)
Compute action matrix-vector products with the FHxc Kernel.
subroutine, public stda_kernel(aop_evects, evects, is_rks_triplets, qs_env, stda_control, stda_env, sub_env, work_matrices)
Compute action matrix-vector products with the sTDA Kernel.
subroutine, public tddfpt_apply_hfx(aop_evects, evects, gs_mos, do_admm, qs_env, work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, work_hmat_asymm, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_apply_energy_diff(aop_evects, evects, s_evects, gs_mos, matrix_ks, tddfpt_control)
Apply orbital energy difference term: Aop_evects(spin,state) += KS(spin) * evects(spin,...
subroutine, public tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, x, res)
...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients – transition cha...
subroutine, public tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, admm_env, hfx_section, x_data, symmetry, recalc_integrals, work_rho_ia_ao, work_hmat, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section)
Write Ritz vectors to a binary restart file.
subroutine, public add_smearing_aterm(qs_env, aop_evects, evects, s_evects, mos_occ, fermia, matrix_s)
...
subroutine, public compute_fermib(qs_env, gs_mos, evals)
...
subroutine, public orthogonalize_smeared_occupation(evects, qs_env, mos_occ, s_c0)
...
pure integer(kind=int_8) function, public tddfpt_total_number_of_states(tddfpt_control, gs_mos)
Compute the number of possible singly excited states (occ -> virt)
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
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...
stores all the informations relevant to an mpi environment
Collection of variables required to evaluate adiabatic TDDFPT kernel.
Type to hold environments for the different kernels.
Parallel (sub)group environment.
Ground state molecular orbitals.
Set of temporary ("work") matrices.