77#include "./base/base_uses.f90"
83 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_tddfpt2_utils'
85 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
87 INTEGER,
PARAMETER,
PRIVATE :: nderivs = 3
88 INTEGER,
PARAMETER,
PRIVATE :: maxspins = 2
108 INTEGER,
INTENT(IN) :: iounit
110 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_init_mos'
112 INTEGER :: handle, ispin, nmo_avail, nmo_occ, &
114 INTEGER,
DIMENSION(2, 2) :: moc, mvt
115 LOGICAL :: print_virtuals_newtonx
116 REAL(kind=
dp),
DIMENSION(:),
POINTER :: evals_virt_spin
120 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
123 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, matrix_s
132 CALL timeset(routinen, handle)
134 CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
135 matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
136 tddfpt_control => dft_control%tddfpt2_control
137 IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. &
138 (tddfpt_control%do_bse_gw_only))
THEN
139 NULLIFY (ks_env, ex_env)
140 CALL get_qs_env(qs_env, exstate_env=ex_env, ks_env=ks_env)
141 CALL dbcsr_copy(matrix_ks(1)%matrix, ex_env%matrix_ks(1)%matrix)
145 cpassert(.NOT.
ASSOCIATED(gs_mos))
147 nspins = dft_control%nspins
148 ALLOCATE (gs_mos(nspins))
152 CALL section_vals_val_get(print_section,
"NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals_newtonx)
157 NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
158 IF (
ASSOCIATED(scf_env))
THEN
159 IF ((scf_env%method ==
ot_method_nr .AND. tddfpt_control%nlumo > 0) .OR. &
160 (scf_env%method ==
ot_method_nr .AND. print_virtuals_newtonx))
THEN
166 CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
167 nmo_virt = min(nmo_virt, nmo_avail - nmo_occ)
170 nmo_virt = tddfpt_control%nlumo - nmo_virt
171 IF (.NOT. print_virtuals_newtonx)
THEN
172 IF (nmo_virt > 0)
THEN
173 ALLOCATE (evals_virt(nspins), mos_virt(nspins))
175 CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
182 IF (
ASSOCIATED(evals_virt))
THEN
183 evals_virt_spin => evals_virt(ispin)%array
185 NULLIFY (evals_virt_spin)
187 IF (
ALLOCATED(mos_virt))
THEN
188 mos_virt_spin => mos_virt(ispin)
190 NULLIFY (mos_virt_spin)
193 nlumo=tddfpt_control%nlumo, &
195 matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
196 mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, &
203 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
204 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
207 WRITE (iounit,
"(T2,A,T36,A)")
"TDDFPT| Molecular Orbitals:", &
208 " Spin AOs Occ Virt Total"
210 WRITE (iounit,
"(T2,A,T37,I4,4I10)")
"TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
211 mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
215 IF (
ASSOCIATED(evals_virt))
THEN
216 DO ispin = 1,
SIZE(evals_virt)
217 IF (
ASSOCIATED(evals_virt(ispin)%array))
DEALLOCATE (evals_virt(ispin)%array)
219 DEALLOCATE (evals_virt)
224 CALL timestop(handle)
250 mos_virt, evals_virt, qs_env)
253 INTEGER,
INTENT(in) :: nlumo
255 INTEGER,
INTENT(in) :: cholesky_method
256 TYPE(
dbcsr_type),
POINTER :: matrix_ks, matrix_s
257 TYPE(
cp_fm_type),
INTENT(IN),
POINTER :: mos_virt
258 REAL(kind=
dp),
DIMENSION(:),
POINTER :: evals_virt
261 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_init_ground_state_mos'
262 REAL(kind=
dp),
PARAMETER :: eps_dp = epsilon(0.0_dp)
264 INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, &
265 irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int
266 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
268 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
269 LOGICAL :: do_eigen, print_phases
270 REAL(kind=
dp) :: element, maxocc
271 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :), &
273 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mo_evals_extended, mo_occ_extended, &
276 ao_mo_virt_fm_struct, wfn_fm_struct
277 TYPE(
cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm, &
279 TYPE(
cp_fm_type),
POINTER :: mo_coeff_extended
285 CALL timeset(routinen, handle)
291 CALL blacs_env%get(para_env=para_env)
293 CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, &
294 nelectron=nelectrons, occupation_numbers=mo_occ_scf)
299 nmo_virt = nao - nmo_occ
301 nmo_virt = min(nmo_virt, nlumo)
304 CALL cp_abort(__location__, &
305 'At least one unoccupied molecular orbital is required to calculate excited states.')
309 IF (
ASSOCIATED(evals_virt))
THEN
310 cpassert(
ASSOCIATED(mos_virt))
311 IF (nmo_virt > nmo_scf - nmo_occ +
SIZE(evals_virt)) do_eigen = .true.
313 IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .true.
317 NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
319 CALL cp_fm_struct_create(ao_mo_occ_fm_struct, template_fmstruct=mo_set%mo_coeff%matrix_struct, &
320 ncol_global=nmo_occ, context=blacs_env)
321 CALL cp_fm_struct_create(ao_mo_virt_fm_struct, template_fmstruct=mo_set%mo_coeff%matrix_struct, &
322 ncol_global=nmo_virt, context=blacs_env)
324 NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
325 ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
327 CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
328 gs_mos%nmo_occ = nmo_occ
330 ALLOCATE (gs_mos%evals_occ(nmo_occ))
331 ALLOCATE (gs_mos%evals_virt(nmo_virt))
332 ALLOCATE (gs_mos%phases_occ(nmo_occ))
333 ALLOCATE (gs_mos%phases_virt(nmo_virt))
336 NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
337 NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
341 CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
342 CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
343 ALLOCATE (mos_extended)
344 CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
345 REAL(nelectrons,
dp), maxocc, flexible_electron_count=0.0_dp)
346 CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name=
"mos-extended")
348 CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
349 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
357 mo_occ_extended(imo) = mo_occ_scf(imo)
359 mo_occ_extended(nmo_scf + 1:) = 0.0_dp
372 cpabort(
'CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
374 cpabort(
'CHOLESKY OFF is not implemented in TDDFT.')
383 cholesky_method_inout = cholesky_method
384 CALL eigensolver(matrix_ks_fm=matrix_ks_fm, mo_set=mos_extended, ortho=ortho_fm, &
385 work=work_fm, cholesky_method=cholesky_method_inout, &
386 do_level_shift=.false., level_shift=0.0_dp, use_jacobi=.false.)
394 CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
395 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
404 CALL cp_fm_to_fm(mo_coeff_extended, work_fm, ncol=nmo_occ, source_start=1, target_start=1)
405 CALL cp_fm_get_info(work_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
406 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
408 ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
409 minrow_neg_array(:) = nao
410 minrow_pos_array(:) = nao
411 sum_sign_array(:) = 0
412 DO icol_local = 1, ncol_local
413 icol_global = col_indices(icol_local)
415 DO irow_local = 1, nrow_local
416 element = my_block(irow_local, icol_local)
419 IF (element >= eps_dp)
THEN
421 ELSE IF (element <= -eps_dp)
THEN
425 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
427 irow_global = row_indices(irow_local)
428 IF (sign_int > 0)
THEN
429 IF (minrow_pos_array(icol_global) > irow_global) &
430 minrow_pos_array(icol_global) = irow_global
431 ELSE IF (sign_int < 0)
THEN
432 IF (minrow_neg_array(icol_global) > irow_global) &
433 minrow_neg_array(icol_global) = irow_global
438 CALL para_env%sum(sum_sign_array)
439 CALL para_env%min(minrow_neg_array)
440 CALL para_env%min(minrow_pos_array)
442 DO icol_local = 1, nmo_occ
443 IF (sum_sign_array(icol_local) > 0)
THEN
445 gs_mos%phases_occ(icol_local) = 1.0_dp
446 ELSE IF (sum_sign_array(icol_local) < 0)
THEN
448 gs_mos%phases_occ(icol_local) = -1.0_dp
451 IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local))
THEN
454 gs_mos%phases_occ(icol_local) = 1.0_dp
457 gs_mos%phases_occ(icol_local) = -1.0_dp
462 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
465 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
466 gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
468 IF (
ASSOCIATED(evals_virt) .AND. (.NOT. do_eigen) .AND. nmo_virt > nmo_scf - nmo_occ)
THEN
469 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_scf - nmo_occ, &
470 source_start=nmo_occ + 1, target_start=1)
471 CALL cp_fm_to_fm(mos_virt, gs_mos%mos_virt, ncol=nmo_virt - (nmo_scf - nmo_occ), &
472 source_start=1, target_start=nmo_scf - nmo_occ + 1)
473 gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
474 gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
476 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
477 gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
480 IF (print_phases)
THEN
486 CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1)
487 CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, &
488 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
490 ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
491 minrow_neg_array(:) = nao
492 minrow_pos_array(:) = nao
493 sum_sign_array(:) = 0
494 DO icol_local = 1, ncol_local
495 icol_global = col_indices(icol_local)
497 DO irow_local = 1, nrow_local
498 element = my_block(irow_local, icol_local)
501 IF (element >= eps_dp)
THEN
503 ELSE IF (element <= -eps_dp)
THEN
507 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
509 irow_global = row_indices(irow_local)
510 IF (sign_int > 0)
THEN
511 IF (minrow_pos_array(icol_global) > irow_global) &
512 minrow_pos_array(icol_global) = irow_global
513 ELSE IF (sign_int < 0)
THEN
514 IF (minrow_neg_array(icol_global) > irow_global) &
515 minrow_neg_array(icol_global) = irow_global
520 CALL para_env%sum(sum_sign_array)
521 CALL para_env%min(minrow_neg_array)
522 CALL para_env%min(minrow_pos_array)
523 DO icol_local = 1, nmo_virt
524 IF (sum_sign_array(icol_local) > 0)
THEN
526 gs_mos%phases_virt(icol_local) = 1.0_dp
527 ELSE IF (sum_sign_array(icol_local) < 0)
THEN
529 gs_mos%phases_virt(icol_local) = -1.0_dp
532 IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local))
THEN
535 gs_mos%phases_virt(icol_local) = 1.0_dp
538 gs_mos%phases_virt(icol_local) = -1.0_dp
542 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
551 DEALLOCATE (mos_extended)
554 CALL timestop(handle)
567 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_release_ground_state_mos'
571 CALL timeset(routinen, handle)
573 IF (
ALLOCATED(gs_mos%phases_occ)) &
574 DEALLOCATE (gs_mos%phases_occ)
576 IF (
ALLOCATED(gs_mos%evals_virt)) &
577 DEALLOCATE (gs_mos%evals_virt)
579 IF (
ALLOCATED(gs_mos%evals_occ)) &
580 DEALLOCATE (gs_mos%evals_occ)
582 IF (
ALLOCATED(gs_mos%phases_virt)) &
583 DEALLOCATE (gs_mos%phases_virt)
585 IF (
ALLOCATED(gs_mos%index_active)) &
586 DEALLOCATE (gs_mos%index_active)
588 IF (
ASSOCIATED(gs_mos%evals_occ_matrix))
THEN
590 DEALLOCATE (gs_mos%evals_occ_matrix)
593 IF (
ASSOCIATED(gs_mos%mos_virt))
THEN
595 DEALLOCATE (gs_mos%mos_virt)
598 IF (
ASSOCIATED(gs_mos%mos_occ))
THEN
600 DEALLOCATE (gs_mos%mos_occ)
603 IF (
ASSOCIATED(gs_mos%mos_active))
THEN
605 DEALLOCATE (gs_mos%mos_active)
608 CALL timestop(handle)
622 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks_oep
624 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_oecorr'
626 INTEGER :: handle, iounit, ispin, nao, nmo_occ, &
631 mo_occ_mo_occ_fm_struct
640 CALL timeset(routinen, handle)
646 CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
647 tddfpt_control => dft_control%tddfpt2_control
651 nspins =
SIZE(gs_mos)
652 NULLIFY (matrix_ks_oep)
653 IF (tddfpt_control%oe_corr /=
oe_none)
THEN
655 WRITE (iounit,
"(1X,A)")
"", &
656 "-------------------------------------------------------------------------------", &
657 "- Orbital Eigenvalue Correction Started -", &
658 "-------------------------------------------------------------------------------"
661 CALL cp_warn(__location__, &
662 "Orbital energy correction potential is an experimental feature. "// &
663 "Use it with extreme care")
668 CALL cp_abort(__location__, &
669 "Implementation of orbital energy correction XC-potentials is "// &
670 "currently incompatible with exact-exchange functionals")
676 CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
682 NULLIFY (xc_fun_empty)
687 IF (dft_control%qs_control%semi_empirical)
THEN
688 cpabort(
"TDDFPT with SE not possible")
689 ELSEIF (dft_control%qs_control%dftb)
THEN
690 cpabort(
"TDDFPT with DFTB not possible")
691 ELSEIF (dft_control%qs_control%xtb)
THEN
692 IF (dft_control%qs_control%xtb_control%do_tblite)
THEN
694 ext_ks_matrix=matrix_ks_oep)
697 ext_ks_matrix=matrix_ks_oep)
701 ext_ks_matrix=matrix_ks_oep)
704 IF (tddfpt_control%oe_corr ==
oe_saop .OR. &
705 tddfpt_control%oe_corr ==
oe_lb .OR. &
706 tddfpt_control%oe_corr ==
oe_gllb)
THEN
708 WRITE (iounit,
"(T2,A)")
" Orbital energy correction of SAOP type "
710 CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
711 ELSE IF (tddfpt_control%oe_corr ==
oe_shift)
THEN
713 WRITE (iounit,
"(T2,A,T71,F10.3)") &
714 " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*
evolt
715 WRITE (iounit,
"(T2,A,T71,F10.3)") &
716 " Open Shell Orbital Eigenvalue Shift [eV] ", tddfpt_control%eos_shift*
evolt
718 CALL ev_shift_operator(qs_env, gs_mos, matrix_ks_oep, &
719 tddfpt_control%ev_shift, tddfpt_control%eos_shift)
721 CALL cp_abort(__location__, &
722 "Unimplemented orbital energy correction potential")
729 NULLIFY (mo_occ_mo_occ_fm_struct)
731 nmo_occ =
SIZE(gs_mos(ispin)%evals_occ)
732 CALL cp_fm_struct_create(mo_occ_mo_occ_fm_struct, nrow_global=nmo_occ, ncol_global=nmo_occ, &
734 ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
735 CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
743 work_fm, ncol=nmo_occ, alpha=1.0_dp, beta=0.0_dp)
744 CALL parallel_gemm(
'T',
'N', nmo_occ, nmo_occ, nao, 1.0_dp, gs_mos(ispin)%mos_occ, work_fm, &
745 0.0_dp, gs_mos(ispin)%evals_occ_matrix)
749 WRITE (iounit,
"(1X,A)") &
750 "-------------------------------------------------------------------------------"
755 CALL timestop(handle)
771 INTEGER(kind=int_8) :: nstates_total
773 INTEGER :: ispin, nspins
776 nspins =
SIZE(gs_mos)
781 nstates_total = nstates_total + &
782 gs_mos(ispin)%nmo_active* &
783 SIZE(gs_mos(ispin)%evals_virt, kind=
int_8)
787 nstates_total = gs_mos(1)%nmo_active* &
788 SIZE(gs_mos(2)%evals_virt, kind=
int_8)
805 SUBROUTINE ev_shift_operator(qs_env, gs_mos, matrix_ks, ev_shift, eos_shift)
811 REAL(kind=
dp),
INTENT(IN) :: ev_shift, eos_shift
813 CHARACTER(len=*),
PARAMETER :: routinen =
'ev_shift_operator'
815 INTEGER :: handle, ispin, n_spins, na, nb, nhomo, &
816 nl, nos, nrow, nu, nvirt
824 CALL timeset(routinen, handle)
826 n_spins =
SIZE(gs_mos)
827 cpassert(n_spins ==
SIZE(matrix_ks))
829 IF (eos_shift /= 0.0_dp .AND. n_spins > 1)
THEN
830 cpabort(
"eos_shift not implemented")
831 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
832 smat => matrix_s(1)%matrix
838 DO ispin = 1, n_spins
839 coeff => gs_mos(ispin)%mos_occ
840 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
841 IF (nhomo == nu)
THEN
850 alpha=-eos_shift, keep_sparsity=.true.)
855 coeff => gs_mos(ispin)%mos_virt
856 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
858 cpassert(nvirt >= nos)
862 alpha=eos_shift, keep_sparsity=.true.)
867 IF (ev_shift /= 0.0_dp)
THEN
868 DO ispin = 1, n_spins
869 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
870 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
871 coeff => gs_mos(ispin)%mos_occ
872 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
876 alpha=-ev_shift, keep_sparsity=.true.)
880 coeff => gs_mos(ispin)%mos_virt
881 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
882 cpassert(nvirt >= nos)
886 alpha=-ev_shift, keep_sparsity=.true.)
893 IF (ev_shift /= 0.0_dp)
THEN
894 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
895 smat => matrix_s(1)%matrix
896 DO ispin = 1, n_spins
897 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
898 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
899 coeff => gs_mos(ispin)%mos_occ
900 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
904 alpha=-ev_shift, keep_sparsity=.true.)
910 IF (eos_shift == 0.0_dp .OR. n_spins == 1)
THEN
911 DO ispin = 1, n_spins
912 IF (
ALLOCATED(gs_mos(ispin)%evals_virt))
THEN
913 gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
923 IF (
ALLOCATED(gs_mos(1)%evals_occ))
THEN
924 gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
926 IF (
ALLOCATED(gs_mos(1)%evals_virt))
THEN
927 gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
929 IF (
ALLOCATED(gs_mos(2)%evals_virt))
THEN
930 gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
931 gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
934 IF (
ALLOCATED(gs_mos(1)%evals_virt))
THEN
935 gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
936 gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
938 IF (
ALLOCATED(gs_mos(2)%evals_occ))
THEN
939 gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
941 IF (
ALLOCATED(gs_mos(2)%evals_virt))
THEN
942 gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
947 CALL timestop(handle)
949 END SUBROUTINE ev_shift_operator
976 fm_pool_ao_mo_active, qs_env, nspins)
977 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(inout) :: evects
978 REAL(kind=
dp),
DIMENSION(:),
INTENT(inout) :: evals
979 INTEGER,
INTENT(in) :: nspins
983 INTEGER,
INTENT(in) :: log_unit
987 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_guess_vectors'
989 CHARACTER(len=5) :: spin_label1, spin_label2
990 INTEGER :: handle, i, imo_occ, imo_virt, ind, ispin, istate, j, jspin, k, no, nstates, &
991 nstates_occ_virt_alpha, nstates_selected, nv, spin1, spin2
992 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: inds
993 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: reverse_index
994 INTEGER,
DIMENSION(maxspins) :: nmo, nmo_occ_avail, nmo_occ_selected, &
996 REAL(kind=
dp) :: e_occ
997 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: e_virt_minus_occ, ev_occ, ev_virt
1000 CALL timeset(routinen, handle)
1002 nstates =
SIZE(evects, 2)
1004 IF (debug_this_module)
THEN
1005 cpassert(nstates > 0)
1006 cpassert(nspins == 1 .OR. nspins == 2)
1012 DO ispin = 1, nspins
1014 nmo_occ_avail(ispin) = gs_mos(ispin)%nmo_active
1015 nmo(ispin) = gs_mos(ispin)%nmo_occ
1018 nmo_occ_selected(ispin) = min(nmo_occ_avail(ispin), nstates)
1019 nmo_virt_selected(ispin) = min(
SIZE(gs_mos(ispin)%evals_virt), nstates)
1025 nstates_selected = dot_product(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
1027 nstates_selected = nmo_occ_selected(1)*nmo_virt_selected(2)
1030 ALLOCATE (inds(nstates_selected))
1031 ALLOCATE (e_virt_minus_occ(nstates_selected))
1036 DO ispin = 1, nspins
1037 no = nmo_occ_selected(ispin)
1038 nv = nmo_virt_selected(ispin)
1039 ALLOCATE (ev_virt(nv), ev_occ(no))
1041 IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. &
1042 (tddfpt_control%do_bse_gw_only))
THEN
1043 ev_virt(1:nv) = ex_env%gw_eigen(nmo(ispin) + 1:nmo(ispin) + nv)
1045 j = nmo_occ_avail(ispin) - i + 1
1046 k = gs_mos(ispin)%index_active(j)
1047 ev_occ(i) = ex_env%gw_eigen(k)
1050 ev_virt(1:nv) = gs_mos(ispin)%evals_virt(1:nv)
1052 j = nmo_occ_avail(ispin) - i + 1
1053 k = gs_mos(ispin)%index_active(j)
1054 ev_occ(i) = gs_mos(ispin)%evals_occ(k)
1058 DO imo_occ = 1, nmo_occ_selected(ispin)
1060 e_occ = ev_occ(imo_occ)
1062 DO imo_virt = 1, nmo_virt_selected(ispin)
1064 e_virt_minus_occ(istate) = ev_virt(imo_virt) - e_occ
1068 DEALLOCATE (ev_virt, ev_occ)
1072 DO imo_occ = 1, nmo_occ_selected(1)
1074 i = gs_mos(1)%nmo_active - imo_occ + 1
1075 k = gs_mos(1)%index_active(i)
1076 e_occ = gs_mos(1)%evals_occ(k)
1078 DO imo_virt = 1, nmo_virt_selected(2)
1080 e_virt_minus_occ(istate) = gs_mos(2)%evals_virt(imo_virt) - e_occ
1085 IF (debug_this_module)
THEN
1086 cpassert(istate == nstates_selected)
1089 CALL sort(e_virt_minus_occ, nstates_selected, inds)
1092 IF (nspins == 1)
THEN
1096 spin_label2 = spin_label1
1101 spin_label1 =
'(alp)'
1102 spin_label2 =
'(bet)'
1107 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
1110 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(2)
1112 IF (log_unit > 0)
THEN
1113 WRITE (log_unit,
"(1X,A)")
"", &
1114 "-------------------------------------------------------------------------------", &
1115 "- TDDFPT Initial Guess -", &
1116 "-------------------------------------------------------------------------------"
1117 WRITE (log_unit,
'(T11,A)')
"State Occupied -> Virtual Excitation"
1118 WRITE (log_unit,
'(T11,A)')
"number orbital orbital energy (eV)"
1119 WRITE (log_unit,
'(1X,79("-"))')
1123 ALLOCATE (reverse_index(i, nspins))
1125 DO ispin = 1, nspins
1126 DO i = 1,
SIZE(gs_mos(ispin)%index_active)
1127 j = gs_mos(ispin)%index_active(i)
1128 reverse_index(j, ispin) = i
1132 DO istate = 1, nstates
1133 IF (
ASSOCIATED(evects(1, istate)%matrix_struct))
THEN
1136 WRITE (log_unit,
'(T7,I8,T28,A19,T60,F14.5)') &
1137 istate,
"*** restarted ***", evals(istate)*
evolt
1142 ind = inds(istate) - 1
1145 IF ((nspins > 1) .AND. (tddfpt_control%spinflip ==
no_sf_tddfpt))
THEN
1146 IF (ind < nstates_occ_virt_alpha)
THEN
1149 spin_label1 =
'(alp)'
1150 spin_label2 =
'(alp)'
1152 ind = ind - nstates_occ_virt_alpha
1155 spin_label1 =
'(bet)'
1156 spin_label2 =
'(bet)'
1162 i = ind/nmo_virt_selected(spin2) + 1
1163 j = nmo_occ_avail(spin1) - i + 1
1164 imo_occ = gs_mos(spin1)%index_active(j)
1165 imo_virt = mod(ind, nmo_virt_selected(spin2)) + 1
1167 evals(istate) = e_virt_minus_occ(istate)
1170 WRITE (log_unit,
'(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
1171 istate, imo_occ, spin_label1, nmo(spin2) + imo_virt, spin_label2, e_virt_minus_occ(istate)*
evolt
1173 DO jspin = 1,
SIZE(evects, 1)
1178 IF (jspin == spin1)
THEN
1181 i = reverse_index(imo_occ, spin1)
1182 CALL cp_fm_to_fm(gs_mos(spin2)%mos_virt, evects(spin1, istate), &
1183 ncol=1, source_start=imo_virt, target_start=i)
1189 DEALLOCATE (reverse_index)
1191 IF (log_unit > 0)
THEN
1192 WRITE (log_unit,
'(/,T7,A,T50,I24)')
'Number of active states:', &
1194 WRITE (log_unit,
"(1X,A)") &
1195 "-------------------------------------------------------------------------------"
1198 DEALLOCATE (e_virt_minus_occ)
1201 CALL timestop(handle)
Handles all functions related to the CELL.
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
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_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
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)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
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
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
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
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_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_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Types for excited states potential energies.
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
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.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix, ext_xc_section)
routine where the real calculations are made: the KS matrix is calculated
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition and initialisation of the mo data type.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name, counter)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
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....
Does all kind of post scf calculations for GPW/GAPW.
subroutine, public make_lumo_gpw(qs_env, scf_env, unoccupied_orbs, unoccupied_evals, nlumo, nlumos)
Gets the LUMOs and their eigenvalues for all spin channels.
module that contains the definitions of the scf types
integer, parameter, public ot_method_nr
subroutine, public tddfpt_release_ground_state_mos(gs_mos)
Release molecular orbitals.
subroutine, public tddfpt_init_ground_state_mos(gs_mos, mo_set, nlumo, blacs_env, cholesky_method, matrix_ks, matrix_s, mos_virt, evals_virt, qs_env)
Generate all virtual molecular orbitals for a given spin by diagonalising the corresponding Kohn-Sham...
subroutine, public tddfpt_oecorr(qs_env, gs_mos, matrix_ks_oep)
Callculate orbital corrected KS matrix for TDDFPT.
subroutine, public tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, fm_pool_ao_mo_active, qs_env, nspins)
Generate missed guess vectors.
subroutine, public tddfpt_init_mos(qs_env, gs_mos, iounit)
Prepare MOs for TDDFPT Calculations.
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)
subroutine, public build_tblite_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
...
All kind of helpful little routines.
Calculate the saop potential.
subroutine, public add_saop_pot(ks_matrix, qs_env, oe_corr)
...
Calculation of KS matrix in xTB Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov JCTC 1...
subroutine, public build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
...
Type defining parameters related to the simulation cell.
represent a pointer to a 1d array
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
to create arrays of pools
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...
Contains information on the excited states energy.
stores all the informations relevant to an mpi environment
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Ground state molecular orbitals.