27 USE dbcsr_api,
ONLY: dbcsr_get_info,&
58 #include "./base/base_uses.f90"
64 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_tddfpt2_eigensolver'
66 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
68 INTEGER,
PARAMETER,
PRIVATE :: nderivs = 3
69 INTEGER,
PARAMETER,
PRIVATE :: maxspins = 2
95 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: evects
96 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(in) :: s_c0_c0t
98 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_orthogonalize_psi1_psi0'
100 INTEGER :: handle, ispin, ivect, nactive, nao, &
102 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct
103 TYPE(cp_fm_type) :: evortho
105 CALL timeset(routinen, handle)
107 nspins =
SIZE(evects, 1)
108 nvects =
SIZE(evects, 2)
112 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
113 nrow_global=nao, ncol_global=nactive)
117 CALL parallel_gemm(
'T',
'N', nao, nactive, nao, 1.0_dp, s_c0_c0t(ispin), &
118 evects(ispin, ivect), 0.0_dp, evortho)
121 CALL cp_fm_release(evortho)
125 CALL timestop(handle)
142 FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm) &
144 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: evects
145 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(in) :: s_c0
146 REAL(kind=
dp),
INTENT(in) :: max_norm
147 LOGICAL :: is_nonortho
149 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_is_nonorthogonal_psi1_psi0'
151 INTEGER :: handle, ispin, ivect, nactive, nao, &
153 REAL(kind=
dp) :: maxabs_val
154 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct, matrix_struct_tmp
155 TYPE(cp_fm_type) :: aortho
157 CALL timeset(routinen, handle)
159 nspins =
SIZE(evects, 1)
160 nvects =
SIZE(evects, 2)
162 is_nonortho = .false.
164 loop:
DO ispin = 1, nspins
166 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
167 nrow_global=nao, ncol_global=nactive)
169 ncol_global=nactive, template_fmstruct=matrix_struct)
174 CALL parallel_gemm(
'T',
'N', nocc, nactive, nao, 1.0_dp, s_c0(ispin), &
175 evects(ispin, ivect), 0.0_dp, aortho)
177 is_nonortho = maxabs_val > max_norm
178 IF (is_nonortho)
THEN
179 CALL cp_fm_release(aortho)
183 CALL cp_fm_release(aortho)
186 CALL timestop(handle)
188 END FUNCTION tddfpt_is_nonorthogonal_psi1_psi0
225 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: evects
226 INTEGER,
INTENT(in) :: nvects_new
227 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: s_evects
228 TYPE(dbcsr_type),
POINTER :: matrix_s
230 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_orthonormalize_psi1_psi1'
232 INTEGER :: handle, ispin, ivect, jvect, nspins, &
233 nvects_old, nvects_total
234 INTEGER,
DIMENSION(maxspins) :: nactive
235 REAL(kind=
dp) :: norm
236 REAL(kind=
dp),
DIMENSION(maxspins) :: weights
238 CALL timeset(routinen, handle)
240 nspins =
SIZE(evects, 1)
241 nvects_total =
SIZE(evects, 2)
242 nvects_old = nvects_total - nvects_new
244 IF (debug_this_module)
THEN
245 cpassert(
SIZE(s_evects, 1) == nspins)
246 cpassert(
SIZE(s_evects, 2) == nvects_total)
247 cpassert(nvects_old >= 0)
251 CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
254 DO jvect = nvects_old + 1, nvects_total
256 DO ivect = 1, jvect - 1
257 CALL cp_fm_trace(evects(:, jvect), s_evects(:, ivect), weights(1:nspins), accurate=.false.)
258 norm = sum(weights(1:nspins))
268 ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
271 CALL cp_fm_trace(evects(:, jvect), s_evects(:, jvect), weights(1:nspins), accurate=.false.)
273 norm = sum(weights(1:nspins))
274 norm = 1.0_dp/sqrt(norm)
282 CALL timestop(handle)
303 SUBROUTINE tddfpt_compute_aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, &
304 matrix_ks, qs_env, kernel_env, &
305 sub_env, work_matrices)
306 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: aop_evects, evects, s_evects
307 TYPE(tddfpt_ground_state_mos),
DIMENSION(:), &
309 TYPE(tddfpt2_control_type),
POINTER :: tddfpt_control
310 TYPE(dbcsr_p_type),
DIMENSION(:),
INTENT(in) :: matrix_ks
311 TYPE(qs_environment_type),
POINTER :: qs_env
312 TYPE(kernel_env_type),
INTENT(in) :: kernel_env
313 TYPE(tddfpt_subgroup_env_type),
INTENT(in) :: sub_env
314 TYPE(tddfpt_work_matrices),
INTENT(inout) :: work_matrices
316 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_Aop_evects'
318 INTEGER :: handle, ispin, ivect, nspins, nvects
319 INTEGER,
DIMENSION(maxspins) :: nmo_occ
320 LOGICAL :: do_admm, do_hfx, do_lri_response, &
321 is_rks_triplets, re_int
322 REAL(kind=
dp) :: rcut, scale
323 TYPE(cp_fm_type) :: fm_dummy
324 TYPE(full_kernel_env_type),
POINTER :: kernel_env_admm_aux
325 TYPE(mp_para_env_type),
POINTER :: para_env
327 CALL timeset(routinen, handle)
329 nspins =
SIZE(evects, 1)
330 nvects =
SIZE(evects, 2)
331 do_hfx = tddfpt_control%do_hfx
332 do_admm = tddfpt_control%do_admm
334 kernel_env_admm_aux => kernel_env%admm_kernel
336 NULLIFY (kernel_env_admm_aux)
338 is_rks_triplets = tddfpt_control%rks_triplets
339 do_lri_response = tddfpt_control%do_lrigpw
341 IF (debug_this_module)
THEN
343 cpassert(
SIZE(aop_evects, 1) == nspins)
344 cpassert(
SIZE(aop_evects, 2) == nvects)
345 cpassert(
SIZE(s_evects, 1) == nspins)
346 cpassert(
SIZE(s_evects, 2) == nvects)
347 cpassert(
SIZE(gs_mos) == nspins)
351 nmo_occ(ispin) =
SIZE(gs_mos(ispin)%evals_occ)
356 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
359 associate(evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect))
360 IF (
ASSOCIATED(evect%matrix_struct))
THEN
361 IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
366 ELSE IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
378 CALL fhxc_kernel(aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, &
379 kernel_env%full_kernel, kernel_env_admm_aux, sub_env, work_matrices, &
380 tddfpt_control%admm_symm, tddfpt_control%admm_xc_correction, &
384 CALL stda_kernel(aop_evects, evects, is_rks_triplets, qs_env, tddfpt_control%stda_control, &
385 kernel_env%stda_kernel, sub_env, work_matrices)
394 cpabort(
"Kernel type not implemented")
397 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
400 associate(aop_evect => aop_evects(ispin, ivect), &
401 work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
402 IF (
ASSOCIATED(aop_evect%matrix_struct))
THEN
403 IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
408 ELSE IF (
ASSOCIATED(work_matrix%matrix_struct))
THEN
419 CALL tddfpt_apply_energy_diff(aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
420 gs_mos=gs_mos, matrix_ks=matrix_ks)
423 IF (tddfpt_control%kernel == tddfpt_kernel_full)
THEN
425 CALL tddfpt_apply_hfx(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
426 qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
427 work_hmat_symm=work_matrices%hfx_hmat_symm, &
428 work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
429 work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
430 work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
431 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda)
THEN
434 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none)
THEN
438 cpabort(
"Kernel type not implemented")
442 IF (tddfpt_control%kernel == tddfpt_kernel_full)
THEN
443 IF (tddfpt_control%do_hfxsr)
THEN
444 re_int = tddfpt_control%hfxsr_re_int
446 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
447 kernel_env%full_kernel%admm_env, &
448 kernel_env%full_kernel%hfxsr_section, &
449 kernel_env%full_kernel%x_data, 1, re_int, &
450 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
451 work_hmat=work_matrices%hfxsr_hmat_symm, &
452 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
454 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
455 kernel_env%full_kernel%admm_env, &
456 kernel_env%full_kernel%hfxsr_section, &
457 kernel_env%full_kernel%x_data, -1, .false., &
458 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
459 work_hmat=work_matrices%hfxsr_hmat_asymm, &
460 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
461 tddfpt_control%hfxsr_re_int = .false.
463 IF (tddfpt_control%do_hfxlr)
THEN
464 rcut = tddfpt_control%hfxlr_rcut
465 scale = tddfpt_control%hfxlr_scale
467 IF (
ALLOCATED(work_matrices%evects_sub))
THEN
468 IF (
ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct))
THEN
469 CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
470 work_matrices%evects_sub(:, ivect), &
471 work_matrices%Aop_evects_sub(:, ivect))
477 CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
478 evects(:, ivect), aop_evects(:, ivect))
486 CALL timestop(handle)
488 END SUBROUTINE tddfpt_compute_aop_evects
505 SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
507 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: ritz_vects, aop_ritz
508 REAL(kind=dp),
DIMENSION(:),
INTENT(out) :: evals
509 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(IN) :: krylov_vects, aop_krylov
510 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: atilde
511 INTEGER,
INTENT(IN) :: nkvo, nkvn
513 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_ritz_vects'
515 INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, &
518 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: at12, at21, at22, evects_atilde
519 TYPE(cp_blacs_env_type),
POINTER :: blacs_env_global
520 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
521 TYPE(cp_fm_type) :: atilde_fm, evects_atilde_fm
523 CALL timeset(routinen, handle)
525 nspins =
SIZE(krylov_vects, 1)
526 nkvs =
SIZE(krylov_vects, 2)
527 nrvs =
SIZE(ritz_vects, 2)
528 cpassert(nkvs == nkvo + nkvn)
530 CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
532 CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
533 CALL cp_fm_create(atilde_fm, fm_struct)
534 CALL cp_fm_create(evects_atilde_fm, fm_struct)
535 CALL cp_fm_struct_release(fm_struct)
538 CALL reallocate(atilde, 1, nkvs, 1, nkvs)
543 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
544 atilde(1:nkvs, 1:nkvs), accurate=.false.)
546 ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
547 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
548 at12, accurate=.false.)
549 atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
550 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
551 at21, accurate=.false.)
552 atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
553 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
554 at22, accurate=.false.)
555 atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
556 DEALLOCATE (at12, at21, at22)
558 atilde = 0.5_dp*(atilde + transpose(atilde))
559 CALL cp_fm_set_submatrix(atilde_fm, atilde)
562 CALL choose_eigv_solver(atilde_fm, evects_atilde_fm, evals(1:nkvs))
564 ALLOCATE (evects_atilde(nkvs, nrvs))
565 CALL cp_fm_get_submatrix(evects_atilde_fm, evects_atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
566 CALL cp_fm_release(evects_atilde_fm)
567 CALL cp_fm_release(atilde_fm)
574 CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
575 CALL cp_fm_set_all(aop_ritz(ispin, irv), 0.0_dp)
579 act = evects_atilde(ikv, irv)
581 CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
582 act, krylov_vects(ispin, ikv))
583 CALL cp_fm_scale_and_add(1.0_dp, aop_ritz(ispin, irv), &
584 act, aop_krylov(ispin, ikv))
590 DEALLOCATE (evects_atilde)
592 CALL timestop(handle)
594 END SUBROUTINE tddfpt_compute_ritz_vects
608 SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
610 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: residual_vects
611 REAL(kind=dp),
DIMENSION(:),
INTENT(in) :: evals
612 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(in) :: ritz_vects, aop_ritz
613 TYPE(tddfpt_ground_state_mos),
DIMENSION(:), &
615 TYPE(dbcsr_type),
POINTER :: matrix_s
617 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_compute_residual_vects'
618 REAL(kind=dp),
PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*epsilon(1.0_dp)
620 INTEGER :: handle, icol_local, irow_local, irv, &
621 ispin, nao, ncols_local, nrows_local, &
623 INTEGER,
DIMENSION(:),
POINTER :: col_indices_local, row_indices_local
624 INTEGER,
DIMENSION(maxspins) :: nactive, nmo_virt
625 REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
626 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :), &
627 POINTER :: weights_ldata
628 TYPE(cp_fm_struct_type),
POINTER :: ao_mo_struct, virt_mo_struct
629 TYPE(cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: awork, vomat
631 CALL timeset(routinen, handle)
633 nspins =
SIZE(residual_vects, 1)
634 nrvs =
SIZE(residual_vects, 2)
637 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
638 ALLOCATE (awork(nspins), vomat(nspins))
640 nmo_virt(ispin) =
SIZE(gs_mos(ispin)%evals_virt)
642 CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
643 ncol_global=nactive(ispin))
644 CALL cp_fm_create(awork(ispin), ao_mo_struct)
646 CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(ispin), &
647 ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
648 CALL cp_fm_create(vomat(ispin), virt_mo_struct)
649 CALL cp_fm_struct_release(virt_mo_struct)
657 CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
658 ncol_local=ncols_local, row_indices=row_indices_local, &
659 col_indices=col_indices_local, local_data=weights_ldata)
662 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
663 ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
664 CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, aop_ritz(ispin, irv))
666 CALL parallel_gemm(
'T',
'N', nmo_virt(ispin), nactive(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
667 awork(ispin), 0.0_dp, vomat(ispin))
669 DO icol_local = 1, ncols_local
670 e_occ_plus_lambda = gs_mos(ispin)%evals_occ(col_indices_local(icol_local)) + lambda
672 DO irow_local = 1, nrows_local
673 eref = gs_mos(ispin)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
677 IF (abs(eref) < threshold) &
678 eref = eref + (1.0_dp - eref_scale)*lambda
680 weights_ldata(irow_local, icol_local) = weights_ldata(irow_local, icol_local)/eref
684 CALL parallel_gemm(
'N',
'N', nao, nactive(ispin), nmo_virt(ispin), 1.0_dp, gs_mos(ispin)%mos_virt, &
685 vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
689 CALL cp_fm_release(awork)
690 CALL cp_fm_release(vomat)
693 CALL timestop(handle)
695 END SUBROUTINE tddfpt_compute_residual_vects
721 matrix_ks, qs_env, kernel_env, &
722 sub_env, logger, iter_unit, energy_unit, &
723 tddfpt_print_section, work_matrices)
RESULT(conv)
724 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(inout) :: evects
725 REAL(kind=dp),
DIMENSION(:),
INTENT(inout) :: evals
726 TYPE(cp_fm_type),
DIMENSION(:, :),
INTENT(inout) :: s_evects
727 TYPE(tddfpt_ground_state_mos),
DIMENSION(:), &
729 TYPE(tddfpt2_control_type),
POINTER :: tddfpt_control
730 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks
731 TYPE(qs_environment_type),
POINTER :: qs_env
732 TYPE(kernel_env_type),
INTENT(in) :: kernel_env
733 TYPE(tddfpt_subgroup_env_type),
INTENT(in) :: sub_env
734 TYPE(cp_logger_type),
POINTER :: logger
735 INTEGER,
INTENT(in) :: iter_unit, energy_unit
736 TYPE(section_vals_type),
POINTER :: tddfpt_print_section
737 TYPE(tddfpt_work_matrices),
INTENT(inout) :: work_matrices
738 REAL(kind=dp) :: conv
740 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_davidson_solver'
742 INTEGER :: handle, ispin, istate, iter, &
743 max_krylov_vects, nspins, nstates, &
744 nstates_conv, nvects_exist, nvects_new
745 INTEGER(kind=int_8) :: nstates_total
746 LOGICAL :: is_nonortho
747 REAL(kind=dp) :: t1, t2
748 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: evals_last
749 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: atilde
750 TYPE(cp_fm_type),
ALLOCATABLE,
DIMENSION(:, :) :: aop_krylov, aop_ritz, krylov_vects, &
752 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
754 CALL timeset(routinen, handle)
756 nspins =
SIZE(gs_mos)
757 nstates = tddfpt_control%nstates
758 nstates_total = tddfpt_total_number_of_states(gs_mos)
760 IF (debug_this_module)
THEN
761 cpassert(
SIZE(evects, 1) == nspins)
762 cpassert(
SIZE(evects, 2) == nstates)
763 cpassert(
SIZE(evals) == nstates)
766 CALL get_qs_env(qs_env, matrix_s=matrix_s)
769 max_krylov_vects = tddfpt_control%nkvs
770 IF (max_krylov_vects < nstates) max_krylov_vects = nstates
771 IF (int(max_krylov_vects, kind=int_8) > nstates_total) max_krylov_vects = int(nstates_total)
773 ALLOCATE (aop_ritz(nspins, nstates))
774 DO istate = 1, nstates
776 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_ritz(ispin, istate))
780 ALLOCATE (evals_last(max_krylov_vects))
781 ALLOCATE (aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), &
782 s_krylov(nspins, max_krylov_vects))
784 DO istate = 1, nstates
786 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate))
787 CALL cp_fm_to_fm(evects(ispin, istate), krylov_vects(ispin, istate))
789 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, s_krylov(ispin, istate))
790 CALL cp_fm_to_fm(s_evects(ispin, istate), s_krylov(ispin, istate))
792 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_krylov(ispin, istate))
801 ALLOCATE (atilde(1, 1))
805 CALL cp_iterate(logger%iter_info, iter_nr_out=iter)
807 CALL tddfpt_compute_aop_evects(aop_evects=aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
808 evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
809 s_evects=s_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
810 gs_mos=gs_mos, tddfpt_control=tddfpt_control, &
811 matrix_ks=matrix_ks, &
812 qs_env=qs_env, kernel_env=kernel_env, &
814 work_matrices=work_matrices)
816 CALL tddfpt_compute_ritz_vects(ritz_vects=evects, aop_ritz=aop_ritz, &
817 evals=evals_last(1:nvects_exist + nvects_new), &
818 krylov_vects=krylov_vects(:, 1:nvects_exist + nvects_new), &
819 aop_krylov=aop_krylov(:, 1:nvects_exist + nvects_new), &
820 atilde=atilde, nkvo=nvects_exist, nkvn=nvects_new)
822 CALL tddfpt_write_restart(evects=evects, evals=evals_last(1:nstates), gs_mos=gs_mos, &
823 logger=logger, tddfpt_print_section=tddfpt_print_section)
825 conv = maxval(abs(evals_last(1:nstates) - evals(1:nstates)))
827 nvects_exist = nvects_exist + nvects_new
828 IF (nvects_exist + nvects_new > max_krylov_vects) &
829 nvects_new = max_krylov_vects - nvects_exist
830 IF (iter >= tddfpt_control%niters) nvects_new = 0
832 IF (conv > tddfpt_control%conv .AND. nvects_new > 0)
THEN
834 DO istate = 1, nvects_new
836 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
837 krylov_vects(ispin, nvects_exist + istate))
838 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
839 s_krylov(ispin, nvects_exist + istate))
840 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
841 aop_krylov(ispin, nvects_exist + istate))
845 CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
846 evals=evals_last(1:nvects_new), &
847 ritz_vects=evects(:, 1:nvects_new), aop_ritz=aop_ritz(:, 1:nvects_new), &
848 gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix)
851 work_matrices%S_C0_C0T)
854 s_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
856 is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
857 work_matrices%S_C0, tddfpt_control%orthogonal_eps)
861 is_nonortho = .false.
865 IF (energy_unit > 0)
THEN
866 WRITE (energy_unit,
'(/,4X,A,T14,A,T36,A)')
"State",
"Exc. energy (eV)",
"Convergence (eV)"
867 DO istate = 1, nstates
868 WRITE (energy_unit,
'(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
869 evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
871 WRITE (energy_unit, *)
872 CALL m_flush(energy_unit)
875 IF (iter_unit > 0)
THEN
877 DO istate = 1, nstates
878 IF (abs(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) &
879 nstates_conv = nstates_conv + 1
882 WRITE (iter_unit,
'(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
883 CALL m_flush(iter_unit)
887 evals(1:nstates) = evals_last(1:nstates)
890 IF (nvects_new == 0 .OR. is_nonortho)
THEN
901 DO istate = nvects_exist + nvects_new, 1, -1
902 DO ispin = nspins, 1, -1
903 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_krylov(ispin, istate))
904 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, s_krylov(ispin, istate))
905 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate))
908 DEALLOCATE (aop_krylov, krylov_vects, s_krylov)
909 DEALLOCATE (evals_last)
911 DO istate = nstates, 1, -1
912 DO ispin = nspins, 1, -1
913 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_ritz(ispin, istate))
916 DEALLOCATE (aop_ritz)
918 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...
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_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 ...
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
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....
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 tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T)
Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
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 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 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)
Compute action matrix-vector products with the FHxc Kernel.
subroutine, public tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks)
Apply orbital energy difference term: Aop_evects(spin,state) += KS(spin) * evects(spin,...
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_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_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_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section)
Write Ritz vectors to a binary restart file.
pure integer(kind=int_8) function, public tddfpt_total_number_of_states(gs_mos)
Compute the number of possible singly excited states (occ -> virt)