75#include "./base/base_uses.f90"
81 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_tddfpt2_utils'
83 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
85 INTEGER,
PARAMETER,
PRIVATE :: nderivs = 3
86 INTEGER,
PARAMETER,
PRIVATE :: maxspins = 2
106 INTEGER,
INTENT(IN) :: iounit
108 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_init_mos'
110 INTEGER :: handle, ispin, nmo_avail, nmo_occ, &
112 INTEGER,
DIMENSION(2, 2) :: moc, mvt
113 LOGICAL :: print_virtuals_newtonx
114 REAL(kind=
dp),
DIMENSION(:),
POINTER :: evals_virt_spin
118 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:), &
121 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, matrix_s
130 CALL timeset(routinen, handle)
132 CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, &
133 matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env)
134 tddfpt_control => dft_control%tddfpt2_control
135 IF (tddfpt_control%do_bse)
THEN
136 NULLIFY (ks_env, ex_env)
137 CALL get_qs_env(qs_env, exstate_env=ex_env, ks_env=ks_env)
138 CALL dbcsr_copy(matrix_ks(1)%matrix, ex_env%matrix_ks(1)%matrix)
142 cpassert(.NOT.
ASSOCIATED(gs_mos))
144 nspins = dft_control%nspins
145 ALLOCATE (gs_mos(nspins))
149 CALL section_vals_val_get(print_section,
"NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals_newtonx)
154 NULLIFY (evals_virt, evals_virt_spin, mos_virt_spin)
155 IF (
ASSOCIATED(scf_env))
THEN
156 IF ((scf_env%method ==
ot_method_nr .AND. tddfpt_control%nlumo > 0) .OR. &
157 (scf_env%method ==
ot_method_nr .AND. print_virtuals_newtonx))
THEN
163 CALL get_mo_set(mos(ispin), nmo=nmo_avail, homo=nmo_occ)
164 nmo_virt = min(nmo_virt, nmo_avail - nmo_occ)
167 nmo_virt = tddfpt_control%nlumo - nmo_virt
168 IF (.NOT. print_virtuals_newtonx)
THEN
169 IF (nmo_virt > 0)
THEN
170 ALLOCATE (evals_virt(nspins), mos_virt(nspins))
172 CALL make_lumo_gpw(qs_env, scf_env, mos_virt, evals_virt, nmo_virt, nmo_avail)
179 IF (
ASSOCIATED(evals_virt))
THEN
180 evals_virt_spin => evals_virt(ispin)%array
182 NULLIFY (evals_virt_spin)
184 IF (
ALLOCATED(mos_virt))
THEN
185 mos_virt_spin => mos_virt(ispin)
187 NULLIFY (mos_virt_spin)
190 nlumo=tddfpt_control%nlumo, &
192 matrix_ks=matrix_ks(ispin)%matrix, matrix_s=matrix_s(1)%matrix, &
193 mos_virt=mos_virt_spin, evals_virt=evals_virt_spin, &
200 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, nrow_global=moc(1, ispin), ncol_global=moc(2, ispin))
201 CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, nrow_global=mvt(1, ispin), ncol_global=mvt(2, ispin))
204 WRITE (iounit,
"(T2,A,T36,A)")
"TDDFPT| Molecular Orbitals:", &
205 " Spin AOs Occ Virt Total"
207 WRITE (iounit,
"(T2,A,T37,I4,4I10)")
"TDDFPT| ", ispin, moc(1, ispin), moc(2, ispin), &
208 mvt(2, ispin), moc(2, ispin) + mvt(2, ispin)
212 IF (
ASSOCIATED(evals_virt))
THEN
213 DO ispin = 1,
SIZE(evals_virt)
214 IF (
ASSOCIATED(evals_virt(ispin)%array))
DEALLOCATE (evals_virt(ispin)%array)
216 DEALLOCATE (evals_virt)
221 CALL timestop(handle)
247 mos_virt, evals_virt, qs_env)
250 INTEGER,
INTENT(in) :: nlumo
252 INTEGER,
INTENT(in) :: cholesky_method
253 TYPE(
dbcsr_type),
POINTER :: matrix_ks, matrix_s
254 TYPE(
cp_fm_type),
INTENT(IN),
POINTER :: mos_virt
255 REAL(kind=
dp),
DIMENSION(:),
POINTER :: evals_virt
258 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_init_ground_state_mos'
259 REAL(kind=
dp),
PARAMETER :: eps_dp = epsilon(0.0_dp)
261 INTEGER :: cholesky_method_inout, handle, icol_global, icol_local, imo, iounit, irow_global, &
262 irow_local, nao, ncol_local, nelectrons, nmo_occ, nmo_scf, nmo_virt, nrow_local, sign_int
263 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
265 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
266 LOGICAL :: do_eigen, print_phases
267 REAL(kind=
dp) :: element, maxocc
268 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :), &
270 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mo_evals_extended, mo_occ_extended, &
273 ao_mo_virt_fm_struct, wfn_fm_struct
274 TYPE(
cp_fm_type) :: matrix_ks_fm, ortho_fm, work_fm, &
276 TYPE(
cp_fm_type),
POINTER :: mo_coeff_extended
282 CALL timeset(routinen, handle)
288 CALL blacs_env%get(para_env=para_env)
290 CALL get_mo_set(mo_set, nao=nao, nmo=nmo_scf, homo=nmo_occ, maxocc=maxocc, &
291 nelectron=nelectrons, occupation_numbers=mo_occ_scf)
296 nmo_virt = nao - nmo_occ
298 nmo_virt = min(nmo_virt, nlumo)
301 CALL cp_abort(__location__, &
302 'At least one unoccupied molecular orbital is required to calculate excited states.')
306 IF (
ASSOCIATED(evals_virt))
THEN
307 cpassert(
ASSOCIATED(mos_virt))
308 IF (nmo_virt > nmo_scf - nmo_occ +
SIZE(evals_virt)) do_eigen = .true.
310 IF (nmo_virt > nmo_scf - nmo_occ) do_eigen = .true.
314 NULLIFY (ao_mo_occ_fm_struct, ao_mo_virt_fm_struct)
315 CALL cp_fm_struct_create(ao_mo_occ_fm_struct, nrow_global=nao, ncol_global=nmo_occ, context=blacs_env)
316 CALL cp_fm_struct_create(ao_mo_virt_fm_struct, nrow_global=nao, ncol_global=nmo_virt, context=blacs_env)
318 NULLIFY (gs_mos%mos_occ, gs_mos%mos_virt, gs_mos%evals_occ_matrix)
319 ALLOCATE (gs_mos%mos_occ, gs_mos%mos_virt)
321 CALL cp_fm_create(gs_mos%mos_virt, ao_mo_virt_fm_struct)
323 ALLOCATE (gs_mos%evals_occ(nmo_occ))
324 ALLOCATE (gs_mos%evals_virt(nmo_virt))
325 ALLOCATE (gs_mos%phases_occ(nmo_occ))
326 ALLOCATE (gs_mos%phases_virt(nmo_virt))
329 NULLIFY (ao_ao_fm_struct, wfn_fm_struct)
330 NULLIFY (mos_extended, mo_coeff_extended, mo_evals_extended, mo_occ_extended)
331 CALL cp_fm_struct_create(ao_ao_fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
335 CALL cp_fm_struct_create(wfn_fm_struct, nrow_global=nao, ncol_global=nmo_occ + nmo_virt, context=blacs_env)
336 ALLOCATE (mos_extended)
337 CALL allocate_mo_set(mos_extended, nao, nmo_occ + nmo_virt, nelectrons, &
338 REAL(nelectrons,
dp), maxocc, flexible_electron_count=0.0_dp)
339 CALL init_mo_set(mos_extended, fm_struct=wfn_fm_struct, name=
"mos-extended")
341 CALL get_mo_set(mos_extended, mo_coeff=mo_coeff_extended, &
342 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
350 mo_occ_extended(imo) = mo_occ_scf(imo)
352 mo_occ_extended(nmo_scf + 1:) = 0.0_dp
365 cpabort(
'CHOLESKY DBCSR_INVERSE is not implemented in TDDFT.')
367 cpabort(
'CHOLESKY OFF is not implemented in TDDFT.')
376 cholesky_method_inout = cholesky_method
377 CALL eigensolver(matrix_ks_fm=matrix_ks_fm, mo_set=mos_extended, ortho=ortho_fm, &
378 work=work_fm, cholesky_method=cholesky_method_inout, &
379 do_level_shift=.false., level_shift=0.0_dp, use_jacobi=.false.)
387 CALL get_mo_set(mo_set, mo_coeff=mo_coeff_extended, &
388 eigenvalues=mo_evals_extended, occupation_numbers=mo_occ_extended)
397 CALL cp_fm_to_fm(mo_coeff_extended, work_fm, ncol=nmo_occ, source_start=1, target_start=1)
398 CALL cp_fm_get_info(work_fm, nrow_local=nrow_local, ncol_local=ncol_local, &
399 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
401 ALLOCATE (minrow_neg_array(nmo_occ), minrow_pos_array(nmo_occ), sum_sign_array(nmo_occ))
402 minrow_neg_array(:) = nao
403 minrow_pos_array(:) = nao
404 sum_sign_array(:) = 0
405 DO icol_local = 1, ncol_local
406 icol_global = col_indices(icol_local)
408 DO irow_local = 1, nrow_local
409 element = my_block(irow_local, icol_local)
412 IF (element >= eps_dp)
THEN
414 ELSE IF (element <= -eps_dp)
THEN
418 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
420 irow_global = row_indices(irow_local)
421 IF (sign_int > 0)
THEN
422 IF (minrow_pos_array(icol_global) > irow_global) &
423 minrow_pos_array(icol_global) = irow_global
424 ELSE IF (sign_int < 0)
THEN
425 IF (minrow_neg_array(icol_global) > irow_global) &
426 minrow_neg_array(icol_global) = irow_global
431 CALL para_env%sum(sum_sign_array)
432 CALL para_env%min(minrow_neg_array)
433 CALL para_env%min(minrow_pos_array)
435 DO icol_local = 1, nmo_occ
436 IF (sum_sign_array(icol_local) > 0)
THEN
438 gs_mos%phases_occ(icol_local) = 1.0_dp
439 ELSE IF (sum_sign_array(icol_local) < 0)
THEN
441 gs_mos%phases_occ(icol_local) = -1.0_dp
444 IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local))
THEN
447 gs_mos%phases_occ(icol_local) = 1.0_dp
450 gs_mos%phases_occ(icol_local) = -1.0_dp
455 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
458 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_occ, ncol=nmo_occ, source_start=1, target_start=1)
459 gs_mos%evals_occ(1:nmo_occ) = mo_evals_extended(1:nmo_occ)
461 IF (
ASSOCIATED(evals_virt) .AND. (.NOT. do_eigen) .AND. nmo_virt > nmo_scf - nmo_occ)
THEN
462 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_scf - nmo_occ, &
463 source_start=nmo_occ + 1, target_start=1)
464 CALL cp_fm_to_fm(mos_virt, gs_mos%mos_virt, ncol=nmo_virt - (nmo_scf - nmo_occ), &
465 source_start=1, target_start=nmo_scf - nmo_occ + 1)
466 gs_mos%evals_virt(1:nmo_scf - nmo_occ) = evals_virt(nmo_occ + 1:nmo_occ + nmo_scf)
467 gs_mos%evals_virt(nmo_scf - nmo_occ + 1:nmo_virt) = evals_virt(1:nmo_virt - (nmo_scf - nmo_occ))
469 CALL cp_fm_to_fm(mo_coeff_extended, gs_mos%mos_virt, ncol=nmo_virt, source_start=nmo_occ + 1, target_start=1)
470 gs_mos%evals_virt(1:nmo_virt) = mo_evals_extended(nmo_occ + 1:nmo_occ + nmo_virt)
473 IF (print_phases)
THEN
479 CALL cp_fm_to_fm(gs_mos%mos_virt, work_fm_virt, ncol=nmo_virt, source_start=1, target_start=1)
480 CALL cp_fm_get_info(work_fm_virt, nrow_local=nrow_local, ncol_local=ncol_local, &
481 row_indices=row_indices, col_indices=col_indices, local_data=my_block)
483 ALLOCATE (minrow_neg_array(nmo_virt), minrow_pos_array(nmo_virt), sum_sign_array(nmo_virt))
484 minrow_neg_array(:) = nao
485 minrow_pos_array(:) = nao
486 sum_sign_array(:) = 0
487 DO icol_local = 1, ncol_local
488 icol_global = col_indices(icol_local)
490 DO irow_local = 1, nrow_local
491 element = my_block(irow_local, icol_local)
494 IF (element >= eps_dp)
THEN
496 ELSE IF (element <= -eps_dp)
THEN
500 sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
502 irow_global = row_indices(irow_local)
503 IF (sign_int > 0)
THEN
504 IF (minrow_pos_array(icol_global) > irow_global) &
505 minrow_pos_array(icol_global) = irow_global
506 ELSE IF (sign_int < 0)
THEN
507 IF (minrow_neg_array(icol_global) > irow_global) &
508 minrow_neg_array(icol_global) = irow_global
513 CALL para_env%sum(sum_sign_array)
514 CALL para_env%min(minrow_neg_array)
515 CALL para_env%min(minrow_pos_array)
516 DO icol_local = 1, nmo_virt
517 IF (sum_sign_array(icol_local) > 0)
THEN
519 gs_mos%phases_virt(icol_local) = 1.0_dp
520 ELSE IF (sum_sign_array(icol_local) < 0)
THEN
522 gs_mos%phases_virt(icol_local) = -1.0_dp
525 IF (minrow_pos_array(icol_local) <= minrow_neg_array(icol_local))
THEN
528 gs_mos%phases_virt(icol_local) = 1.0_dp
531 gs_mos%phases_virt(icol_local) = -1.0_dp
535 DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
544 DEALLOCATE (mos_extended)
547 CALL timestop(handle)
560 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_release_ground_state_mos'
564 CALL timeset(routinen, handle)
566 IF (
ALLOCATED(gs_mos%phases_occ)) &
567 DEALLOCATE (gs_mos%phases_occ)
569 IF (
ALLOCATED(gs_mos%evals_virt)) &
570 DEALLOCATE (gs_mos%evals_virt)
572 IF (
ALLOCATED(gs_mos%evals_occ)) &
573 DEALLOCATE (gs_mos%evals_occ)
575 IF (
ALLOCATED(gs_mos%phases_virt)) &
576 DEALLOCATE (gs_mos%phases_virt)
578 IF (
ASSOCIATED(gs_mos%evals_occ_matrix))
THEN
580 DEALLOCATE (gs_mos%evals_occ_matrix)
583 IF (
ASSOCIATED(gs_mos%mos_virt))
THEN
585 DEALLOCATE (gs_mos%mos_virt)
588 IF (
ASSOCIATED(gs_mos%mos_occ))
THEN
590 DEALLOCATE (gs_mos%mos_occ)
593 CALL timestop(handle)
606 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks_oep
608 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_oecorr'
610 INTEGER :: handle, iounit, ispin, nao, nmo_occ, &
615 mo_occ_mo_occ_fm_struct
624 CALL timeset(routinen, handle)
630 CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, matrix_ks=matrix_ks)
631 tddfpt_control => dft_control%tddfpt2_control
635 nspins =
SIZE(gs_mos)
636 NULLIFY (matrix_ks_oep)
637 IF (tddfpt_control%oe_corr /=
oe_none)
THEN
639 WRITE (iounit,
"(1X,A)")
"", &
640 "-------------------------------------------------------------------------------", &
641 "- Orbital Eigenvalue Correction Started -", &
642 "-------------------------------------------------------------------------------"
645 CALL cp_warn(__location__, &
646 "Orbital energy correction potential is an experimental feature. "// &
647 "Use it with extreme care")
652 CALL cp_abort(__location__, &
653 "Implementation of orbital energy correction XC-potentials is "// &
654 "currently incompatible with exact-exchange functionals")
660 CALL dbcsr_copy(matrix_ks_oep(ispin)%matrix, matrix_ks(ispin)%matrix)
666 NULLIFY (xc_fun_empty)
671 IF (dft_control%qs_control%semi_empirical)
THEN
672 cpabort(
"TDDFPT with SE not possible")
673 ELSEIF (dft_control%qs_control%dftb)
THEN
674 cpabort(
"TDDFPT with DFTB not possible")
675 ELSEIF (dft_control%qs_control%xtb)
THEN
677 ext_ks_matrix=matrix_ks_oep)
680 ext_ks_matrix=matrix_ks_oep)
683 IF (tddfpt_control%oe_corr ==
oe_saop .OR. &
684 tddfpt_control%oe_corr ==
oe_lb .OR. &
685 tddfpt_control%oe_corr ==
oe_gllb)
THEN
687 WRITE (iounit,
"(T2,A)")
" Orbital energy correction of SAOP type "
689 CALL add_saop_pot(matrix_ks_oep, qs_env, tddfpt_control%oe_corr)
690 ELSE IF (tddfpt_control%oe_corr ==
oe_shift)
THEN
692 WRITE (iounit,
"(T2,A,T71,F10.3)") &
693 " Virtual Orbital Eigenvalue Shift [eV] ", tddfpt_control%ev_shift*
evolt
694 WRITE (iounit,
"(T2,A,T71,F10.3)") &
695 " Open Shell Orbital Eigenvalue Shift [eV] ", tddfpt_control%eos_shift*
evolt
697 CALL ev_shift_operator(qs_env, gs_mos, matrix_ks_oep, &
698 tddfpt_control%ev_shift, tddfpt_control%eos_shift)
700 CALL cp_abort(__location__, &
701 "Unimplemented orbital energy correction potential")
708 NULLIFY (mo_occ_mo_occ_fm_struct)
710 nmo_occ =
SIZE(gs_mos(ispin)%evals_occ)
711 CALL cp_fm_struct_create(mo_occ_mo_occ_fm_struct, nrow_global=nmo_occ, ncol_global=nmo_occ, &
713 ALLOCATE (gs_mos(ispin)%evals_occ_matrix)
714 CALL cp_fm_create(gs_mos(ispin)%evals_occ_matrix, mo_occ_mo_occ_fm_struct)
722 work_fm, ncol=nmo_occ, alpha=1.0_dp, beta=0.0_dp)
723 CALL parallel_gemm(
'T',
'N', nmo_occ, nmo_occ, nao, 1.0_dp, gs_mos(ispin)%mos_occ, work_fm, &
724 0.0_dp, gs_mos(ispin)%evals_occ_matrix)
728 WRITE (iounit,
"(1X,A)") &
729 "-------------------------------------------------------------------------------"
734 CALL timestop(handle)
750 INTEGER(kind=int_8) :: nstates_total
752 INTEGER :: ispin, nspins
755 nspins =
SIZE(gs_mos)
760 nstates_total = nstates_total + &
761 SIZE(gs_mos(ispin)%evals_occ, kind=
int_8)* &
762 SIZE(gs_mos(ispin)%evals_virt, kind=
int_8)
766 nstates_total =
SIZE(gs_mos(1)%evals_occ, kind=
int_8)* &
767 SIZE(gs_mos(2)%evals_virt, kind=
int_8)
784 SUBROUTINE ev_shift_operator(qs_env, gs_mos, matrix_ks, ev_shift, eos_shift)
790 REAL(kind=
dp),
INTENT(IN) :: ev_shift, eos_shift
792 CHARACTER(len=*),
PARAMETER :: routinen =
'ev_shift_operator'
794 INTEGER :: handle, ispin, n_spins, na, nb, nhomo, &
795 nl, nos, nrow, nu, nvirt
803 CALL timeset(routinen, handle)
805 n_spins =
SIZE(gs_mos)
806 cpassert(n_spins ==
SIZE(matrix_ks))
808 IF (eos_shift /= 0.0_dp .AND. n_spins > 1)
THEN
809 cpabort(
"eos_shift not implemented")
810 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
811 smat => matrix_s(1)%matrix
817 DO ispin = 1, n_spins
818 coeff => gs_mos(ispin)%mos_occ
819 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
820 IF (nhomo == nu)
THEN
829 alpha=-eos_shift, keep_sparsity=.true.)
834 coeff => gs_mos(ispin)%mos_virt
835 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
837 cpassert(nvirt >= nos)
841 alpha=eos_shift, keep_sparsity=.true.)
846 IF (ev_shift /= 0.0_dp)
THEN
847 DO ispin = 1, n_spins
848 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
849 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
850 coeff => gs_mos(ispin)%mos_occ
851 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
855 alpha=-ev_shift, keep_sparsity=.true.)
859 coeff => gs_mos(ispin)%mos_virt
860 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nvirt)
861 cpassert(nvirt >= nos)
865 alpha=-ev_shift, keep_sparsity=.true.)
872 IF (ev_shift /= 0.0_dp)
THEN
873 CALL get_qs_env(qs_env, mos=mos, matrix_s=matrix_s)
874 smat => matrix_s(1)%matrix
875 DO ispin = 1, n_spins
876 CALL dbcsr_add(matrix_ks(ispin)%matrix, smat, &
877 alpha_scalar=1.0_dp, beta_scalar=ev_shift)
878 coeff => gs_mos(ispin)%mos_occ
879 CALL cp_fm_get_info(coeff, matrix_struct=fmstruct, ncol_global=nhomo)
883 alpha=-ev_shift, keep_sparsity=.true.)
889 IF (eos_shift == 0.0_dp .OR. n_spins == 1)
THEN
890 DO ispin = 1, n_spins
891 IF (
ALLOCATED(gs_mos(ispin)%evals_virt))
THEN
892 gs_mos(ispin)%evals_virt(:) = gs_mos(ispin)%evals_virt(:) + ev_shift
902 IF (
ALLOCATED(gs_mos(1)%evals_occ))
THEN
903 gs_mos(1)%evals_occ(nl + 1:nu) = gs_mos(1)%evals_occ(nl + 1:nu) - eos_shift
905 IF (
ALLOCATED(gs_mos(1)%evals_virt))
THEN
906 gs_mos(1)%evals_virt(:) = gs_mos(1)%evals_virt(:) + ev_shift
908 IF (
ALLOCATED(gs_mos(2)%evals_virt))
THEN
909 gs_mos(2)%evals_virt(1:nos) = gs_mos(2)%evals_virt(1:nos) + eos_shift
910 gs_mos(2)%evals_virt(nos + 1:) = gs_mos(2)%evals_virt(nos + 1:) + ev_shift
913 IF (
ALLOCATED(gs_mos(1)%evals_virt))
THEN
914 gs_mos(1)%evals_virt(1:nos) = gs_mos(1)%evals_virt(1:nos) + eos_shift
915 gs_mos(1)%evals_virt(nos + 1:) = gs_mos(1)%evals_virt(nos + 1:) + ev_shift
917 IF (
ALLOCATED(gs_mos(2)%evals_occ))
THEN
918 gs_mos(2)%evals_occ(nl + 1:nu) = gs_mos(2)%evals_occ(nl + 1:nu) - eos_shift
920 IF (
ALLOCATED(gs_mos(2)%evals_virt))
THEN
921 gs_mos(2)%evals_virt(:) = gs_mos(2)%evals_virt(:) + ev_shift
926 CALL timestop(handle)
928 END SUBROUTINE ev_shift_operator
954 TYPE(
cp_fm_type),
DIMENSION(:, :),
INTENT(inout) :: evects
955 REAL(kind=
dp),
DIMENSION(:),
INTENT(inout) :: evals
956 INTEGER,
INTENT(in) :: nspins
959 INTEGER,
INTENT(in) :: log_unit
963 CHARACTER(LEN=*),
PARAMETER :: routinen =
'tddfpt_guess_vectors'
965 CHARACTER(len=5) :: spin_label1, spin_label2
966 INTEGER :: handle, imo_occ, imo_virt, ind, ispin, istate, jspin, nstates, &
967 nstates_occ_virt_alpha, nstates_selected, spin1, spin2, total
968 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: inds
969 INTEGER,
DIMENSION(maxspins) :: nmo_occ_avail, nmo_occ_selected, &
971 REAL(kind=
dp) :: e_occ
972 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: e_virt_minus_occ, gw_occ, gw_virt
976 CALL timeset(routinen, handle)
978 nstates =
SIZE(evects, 2)
980 IF (debug_this_module)
THEN
981 cpassert(nstates > 0)
982 cpassert(nspins == 1 .OR. nspins == 2)
990 nmo_occ_avail(ispin) =
SIZE(gs_mos(ispin)%evals_occ)
993 nmo_occ_selected(ispin) = min(nmo_occ_avail(ispin), nstates)
994 nmo_virt_selected(ispin) = min(
SIZE(gs_mos(ispin)%evals_virt), nstates)
996 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct_evects(ispin)%struct)
1002 nstates_selected = dot_product(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
1004 nstates_selected = nmo_occ_selected(1)*nmo_virt_selected(2)
1007 ALLOCATE (inds(nstates_selected))
1008 ALLOCATE (e_virt_minus_occ(nstates_selected))
1013 DO ispin = 1, nspins
1015 IF (tddfpt_control%do_bse)
THEN
1016 ALLOCATE (gw_virt(nmo_virt_selected(ispin)), gw_occ(nmo_occ_selected(ispin)))
1017 total = nmo_virt_selected(ispin) + nmo_occ_selected(ispin)
1018 gw_virt(1:nmo_virt_selected(ispin)) = ex_env%gw_eigen((nmo_occ_selected(ispin) + 1):total)
1019 gw_occ(1:nmo_occ_selected(ispin)) = ex_env%gw_eigen(1:nmo_occ_selected(ispin))
1022 DO imo_occ = 1, nmo_occ_selected(ispin)
1024 IF (.NOT. tddfpt_control%do_bse)
THEN
1025 e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1)
1027 e_occ = gw_occ(nmo_occ_avail(ispin) - imo_occ + 1)
1030 DO imo_virt = 1, nmo_virt_selected(ispin)
1032 IF (.NOT. tddfpt_control%do_bse)
THEN
1033 e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ
1035 e_virt_minus_occ(istate) = gw_virt(imo_virt) - e_occ
1040 IF (tddfpt_control%do_bse)
THEN
1041 DEALLOCATE (gw_virt, gw_occ)
1046 DO imo_occ = 1, nmo_occ_selected(1)
1048 e_occ = gs_mos(1)%evals_occ(nmo_occ_avail(1) - imo_occ + 1)
1050 DO imo_virt = 1, nmo_virt_selected(2)
1052 e_virt_minus_occ(istate) = gs_mos(2)%evals_virt(imo_virt) - e_occ
1057 IF (debug_this_module)
THEN
1058 cpassert(istate == nstates_selected)
1061 CALL sort(e_virt_minus_occ, nstates_selected, inds)
1064 IF (nspins == 1)
THEN
1068 spin_label2 = spin_label1
1070 ELSE IF (tddfpt_control%spinflip .NE.
no_sf_tddfpt)
THEN
1073 spin_label1 =
'(alp)'
1074 spin_label2 =
'(bet)'
1079 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
1082 nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(2)
1084 IF (log_unit > 0)
THEN
1085 WRITE (log_unit,
"(1X,A)")
"", &
1086 "-------------------------------------------------------------------------------", &
1087 "- TDDFPT Initial Guess -", &
1088 "-------------------------------------------------------------------------------"
1089 WRITE (log_unit,
'(T11,A)')
"State Occupied -> Virtual Excitation"
1090 WRITE (log_unit,
'(T11,A)')
"number orbital orbital energy (eV)"
1091 WRITE (log_unit,
'(1X,79("-"))')
1094 DO istate = 1, nstates
1096 IF (
ASSOCIATED(evects(1, istate)%matrix_struct))
THEN
1098 WRITE (log_unit,
'(T7,I8,T28,A19,T60,F14.5)') &
1099 istate,
"*** restarted ***", evals(istate)*
evolt
1103 ind = inds(istate) - 1
1106 IF ((nspins > 1) .AND. (tddfpt_control%spinflip ==
no_sf_tddfpt))
THEN
1107 IF (ind < nstates_occ_virt_alpha)
THEN
1110 spin_label1 =
'(alp)'
1111 spin_label2 =
'(alp)'
1113 ind = ind - nstates_occ_virt_alpha
1116 spin_label1 =
'(bet)'
1117 spin_label2 =
'(bet)'
1123 imo_occ = nmo_occ_avail(spin1) - ind/nmo_virt_selected(spin2)
1124 imo_virt = mod(ind, nmo_virt_selected(spin2)) + 1
1126 evals(istate) = e_virt_minus_occ(istate)
1129 WRITE (log_unit,
'(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
1130 istate, imo_occ, spin_label1, nmo_occ_avail(spin2) + imo_virt, spin_label2, e_virt_minus_occ(istate)*
evolt
1132 DO jspin = 1,
SIZE(evects, 1)
1134 CALL cp_fm_create(evects(jspin, istate), fm_struct_evects(jspin)%struct)
1137 IF (jspin == spin1) &
1140 CALL cp_fm_to_fm(gs_mos(spin2)%mos_virt, evects(spin1, istate), &
1141 ncol=1, source_start=imo_virt, target_start=imo_occ)
1146 IF (log_unit > 0)
THEN
1148 WRITE (log_unit,
"(1X,A)") &
1149 "-------------------------------------------------------------------------------"
1152 DEALLOCATE (e_virt_minus_occ)
1155 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,...
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_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_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 ...
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, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
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)
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, 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, 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 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.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
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 eigenvalues for the lumos.
module that contains the definitions of the scf types
integer, parameter, public ot_method_nr
subroutine, public tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, qs_env, nspins)
Generate missed guess vectors.
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_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)
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...
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.