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")
465 rcut = tddfpt_control%hfxlr_rcut
466 scale = tddfpt_control%hfxlr_scale
468 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
469 IF (
ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct))
THEN
471 work_matrices%evects_sub(:, ivect), &
472 work_matrices%Aop_evects_sub(:, ivect))
476 evects(:, ivect), aop_evects(:, ivect))
481 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
483 DO ispin = 1,
SIZE(evects, 1)
484 associate(aop_evect => aop_evects(ispin, ivect), &
485 work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
486 IF (
ASSOCIATED(aop_evect%matrix_struct))
THEN
487 IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
492 ELSE IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
503 IF ((.NOT. do_bse) .AND. (.NOT. do_bse_gw_only))
THEN
504 CALL tddfpt_apply_energy_diff(aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
505 gs_mos=gs_mos, matrix_ks=matrix_ks, tddfpt_control=tddfpt_control)
507 CALL zeroth_order_gw(qs_env=qs_env, aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
508 gs_mos=gs_mos, matrix_s=matrix_s, matrix_ks=matrix_ks)
512 IF (tddfpt_control%do_smearing)
THEN
515 CALL add_smearing_aterm(qs_env, aop_evects(ispin, ivect), evects(ispin, ivect), &
516 s_evects(ispin, ivect), gs_mos(ispin)%mos_occ, &
517 tddfpt_control%smeared_occup(ispin)%fermia, matrix_s)
523 IF (tddfpt_control%kernel == tddfpt_kernel_full)
THEN
525 CALL tddfpt_apply_hfx(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
526 qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
527 work_hmat_symm=work_matrices%hfx_hmat_symm, &
528 work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
529 work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
530 work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
531 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda)
THEN
534 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none)
THEN
538 cpabort(
"Kernel type not implemented")
542 IF (tddfpt_control%kernel == tddfpt_kernel_full)
THEN
543 IF (tddfpt_control%do_hfxsr)
THEN
544 re_int = tddfpt_control%hfxsr_re_int
546 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
547 kernel_env%full_kernel%admm_env, &
548 kernel_env%full_kernel%hfxsr_section, &
549 kernel_env%full_kernel%x_data, 1, re_int, &
550 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
551 work_hmat=work_matrices%hfxsr_hmat_symm, &
552 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
554 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
555 kernel_env%full_kernel%admm_env, &
556 kernel_env%full_kernel%hfxsr_section, &
557 kernel_env%full_kernel%x_data, -1, .false., &
558 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
559 work_hmat=work_matrices%hfxsr_hmat_asymm, &
560 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
561 tddfpt_control%hfxsr_re_int = .false.
564 IF ((do_bse) .OR. (do_bse_w_only))
THEN
567 CALL tddfpt_apply_bse(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env, &
573 CALL timestop(handle)
575 END SUBROUTINE tddfpt_compute_aop_evects
592 SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
594 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: ritz_vects, aop_ritz
595 REAL(kind=dp),
DIMENSION(:),
INTENT(out) :: evals
596 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: krylov_vects, aop_krylov
597 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: atilde
598 INTEGER,
INTENT(IN) :: nkvo, nkvn
600 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_ritz_vects'
602 INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, &
605 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: at12, at21, at22, evects_atilde
606 TYPE(cp_blacs_env_type),
POINTER :: blacs_env_global
607 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
608 TYPE(cp_fm_type) :: atilde_fm, evects_atilde_fm
610 CALL timeset(routinen, handle)
612 nspins =
SIZE(krylov_vects, 1)
613 nkvs =
SIZE(krylov_vects, 2)
614 nrvs =
SIZE(ritz_vects, 2)
615 cpassert(nkvs == nkvo + nkvn)
617 CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
619 CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
620 CALL cp_fm_create(atilde_fm, fm_struct, set_zero=.true.)
621 CALL cp_fm_create(evects_atilde_fm, fm_struct, set_zero=.true.)
622 CALL cp_fm_struct_release(fm_struct)
625 CALL reallocate(atilde, 1, nkvs, 1, nkvs)
630 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
631 atilde(1:nkvs, 1:nkvs), accurate=.false.)
633 ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
634 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
635 at12, accurate=.false.)
636 atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
637 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
638 at21, accurate=.false.)
639 atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
640 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
641 at22, accurate=.false.)
642 atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
643 DEALLOCATE (at12, at21, at22)
645 atilde = 0.5_dp*(atilde + transpose(atilde))
646 CALL cp_fm_set_submatrix(atilde_fm, atilde)
649 CALL choose_eigv_solver(atilde_fm, evects_atilde_fm, evals(1:nkvs))
651 ALLOCATE (evects_atilde(nkvs, nrvs))
652 CALL cp_fm_get_submatrix(evects_atilde_fm, evects_atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
653 CALL cp_fm_release(evects_atilde_fm)
654 CALL cp_fm_release(atilde_fm)
661 CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
662 CALL cp_fm_set_all(aop_ritz(ispin, irv), 0.0_dp)
666 act = evects_atilde(ikv, irv)
668 CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
669 act, krylov_vects(ispin, ikv))
670 CALL cp_fm_scale_and_add(1.0_dp, aop_ritz(ispin, irv), &
671 act, aop_krylov(ispin, ikv))
677 DEALLOCATE (evects_atilde)
679 CALL timestop(handle)
681 END SUBROUTINE tddfpt_compute_ritz_vects
696 SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
697 matrix_s, tddfpt_control)
698 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: residual_vects
699 REAL(kind=dp),
DIMENSION(:),
INTENT(in) :: evals
700 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: ritz_vects, aop_ritz
701 TYPE(tddfpt_ground_state_mos),
DIMENSION(:), &
703 TYPE(dbcsr_type),
POINTER :: matrix_s
704 TYPE(tddfpt2_control_type),
POINTER :: tddfpt_control
706 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_residual_vects'
707 REAL(kind=dp),
PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*epsilon(1.0_dp)
709 INTEGER :: handle, ica, icb, icol_local, &
710 irow_local, irv, ispin, nao, &
711 ncols_local, nrows_local, nrvs, &
712 nspins, spin2, spinflip
713 INTEGER,
DIMENSION(:),
POINTER :: col_indices_local, row_indices_local
714 INTEGER,
DIMENSION(maxspins) :: nactive, nmo_virt
715 REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
716 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :), &
717 POINTER :: weights_ldata
718 TYPE(cp_fm_struct_type),
POINTER :: ao_mo_struct, virt_mo_struct
719 TYPE(cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: awork, vomat
721 CALL timeset(routinen, handle)
723 nspins =
SIZE(residual_vects, 1)
724 nrvs =
SIZE(residual_vects, 2)
725 spinflip = tddfpt_control%spinflip
728 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
729 ALLOCATE (awork(nspins), vomat(nspins))
731 IF (spinflip == no_sf_tddfpt)
THEN
736 nmo_virt(spin2) =
SIZE(gs_mos(spin2)%evals_virt)
738 CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
739 ncol_global=nactive(ispin))
740 CALL cp_fm_create(awork(ispin), ao_mo_struct)
742 CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(spin2), &
743 ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
744 CALL cp_fm_create(vomat(ispin), virt_mo_struct)
745 CALL cp_fm_struct_release(virt_mo_struct)
753 IF (spinflip == no_sf_tddfpt)
THEN
758 CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
759 ncol_local=ncols_local, row_indices=row_indices_local, &
760 col_indices=col_indices_local, local_data=weights_ldata)
763 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
764 ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
765 CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, aop_ritz(ispin, irv))
767 CALL parallel_gemm(
'T',
'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, gs_mos(spin2)%mos_virt, &
768 awork(ispin), 0.0_dp, vomat(ispin))
771 DO icol_local = 1, ncols_local
772 ica = col_indices_local(icol_local)
773 icb = gs_mos(ispin)%index_active(ica)
774 e_occ_plus_lambda = gs_mos(ispin)%evals_occ(icb) + lambda
776 DO irow_local = 1, nrows_local
777 eref = gs_mos(spin2)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
781 IF (abs(eref) < threshold)
THEN
782 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)
THEN
934 nvects_new = max_krylov_vects - nvects_exist
936 IF (iter >= tddfpt_control%niters) nvects_new = 0
938 IF (conv > tddfpt_control%conv .AND. nvects_new > 0)
THEN
940 DO istate = 1, nvects_new
942 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
943 krylov_vects(ispin, nvects_exist + istate))
944 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
945 s_krylov(ispin, nvects_exist + istate))
946 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
947 aop_krylov(ispin, nvects_exist + istate))
951 CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
952 evals=evals_last(1:nvects_new), &
953 ritz_vects=evects(:, 1:nvects_new), aop_ritz=aop_ritz(:, 1:nvects_new), &
954 gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix, tddfpt_control=tddfpt_control)
957 work_matrices%S_C0_C0T, qs_env, &
958 gs_mos, evals(1:nstates), tddfpt_control, work_matrices%S_C0)
961 s_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
963 is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
964 work_matrices%S_C0, tddfpt_control%orthogonal_eps, &
965 tddfpt_control%spinflip)
969 is_nonortho = .false.
973 IF (energy_unit > 0)
THEN
974 WRITE (energy_unit,
'(/,4X,A,T14,A,T36,A)')
"State",
"Exc. energy (eV)",
"Convergence (eV)"
975 DO istate = 1, nstates
976 WRITE (energy_unit,
'(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
977 evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
979 WRITE (energy_unit, *)
980 CALL m_flush(energy_unit)
983 IF (iter_unit > 0)
THEN
985 DO istate = 1, nstates
986 IF (abs(evals_last(istate) - evals(istate)) <= tddfpt_control%conv)
THEN
987 nstates_conv = nstates_conv + 1
991 WRITE (iter_unit,
'(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
992 CALL m_flush(iter_unit)
996 evals(1:nstates) = evals_last(1:nstates)
999 IF (nvects_new == 0 .OR. is_nonortho)
THEN
1003 evals(1:nstates), tddfpt_control, work_matrices%S_C0)
1012 DO istate = nvects_exist + nvects_new, 1, -1
1013 DO ispin = nspins, 1, -1
1014 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_krylov(ispin, istate))
1015 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, s_krylov(ispin, istate))
1016 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
1019 DEALLOCATE (aop_krylov, krylov_vects, s_krylov)
1020 DEALLOCATE (evals_last)
1022 DO istate = nstates, 1, -1
1023 DO ispin = nspins, 1, -1
1024 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_ritz(ispin, istate))
1027 DEALLOCATE (aop_ritz)
1029 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_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_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,...
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, 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.
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.