28 USE dbcsr_api,
ONLY: dbcsr_add,&
31 dbcsr_deallocate_matrix,&
55 lri_environment_type,&
78 integrate_v_rspace_diagonal,&
79 integrate_v_rspace_one_center
107 #include "./base/base_uses.f90"
119 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_linres_kernel'
134 TYPE(qs_environment_type),
POINTER :: qs_env
135 TYPE(qs_p_env_type) :: p_env
136 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(IN) :: c0, av
138 INTEGER :: ispin, ncol
139 TYPE(dft_control_type),
POINTER :: dft_control
141 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
142 IF (dft_control%qs_control%semi_empirical)
THEN
143 cpabort(
"Linear response not available with SE methods")
144 ELSEIF (dft_control%qs_control%dftb)
THEN
145 cpabort(
"Linear response not available with DFTB")
146 ELSEIF (dft_control%qs_control%xtb)
THEN
147 CALL apply_op_2_xtb(qs_env, p_env)
149 CALL apply_op_2_dft(qs_env, p_env)
155 DO ispin = 1,
SIZE(c0)
160 ncol=ncol, alpha=1.0_dp, beta=1.0_dp)
170 SUBROUTINE apply_op_2_dft(qs_env, p_env)
171 TYPE(qs_environment_type),
POINTER :: qs_env
172 TYPE(qs_p_env_type) :: p_env
174 CHARACTER(len=*),
PARAMETER :: routinen =
'apply_op_2_dft'
176 INTEGER :: handle, ikind, ispin, nkind, ns, nspins
177 LOGICAL :: deriv2_analytic, gapw, gapw_xc, &
179 REAL(kind=
dp) :: alpha, ekin_mol, energy_hartree, &
181 TYPE(admm_type),
POINTER :: admm_env
182 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
183 TYPE(cp_logger_type),
POINTER :: logger
184 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: k1mat, matrix_s, rho1_ao, rho_ao
185 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, psmat
186 TYPE(dft_control_type),
POINTER :: dft_control
187 TYPE(kg_environment_type),
POINTER :: kg_env
188 TYPE(linres_control_type),
POINTER :: linres_control
189 TYPE(lri_density_type),
POINTER :: lri_density
190 TYPE(lri_environment_type),
POINTER :: lri_env
191 TYPE(lri_kind_type),
DIMENSION(:),
POINTER :: lri_v_int
192 TYPE(mp_para_env_type),
POINTER :: para_env
193 TYPE(pw_c1d_gs_type) :: rho1_tot_gspace, v_hartree_gspace
194 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho1_g
195 TYPE(pw_env_type),
POINTER :: pw_env
196 TYPE(pw_poisson_type),
POINTER :: poisson_env
197 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
198 TYPE(pw_r3d_rs_type) :: v_hartree_rspace
199 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho1_r, rho_r, tau1_r, v_rspace_new, &
201 TYPE(qs_kpp1_env_type),
POINTER :: kpp1_env
202 TYPE(qs_ks_env_type),
POINTER :: ks_env
203 TYPE(qs_rho_type),
POINTER :: rho, rho0, rho1, rho1_xc, rho1a, &
205 TYPE(rho_atom_type),
DIMENSION(:),
POINTER :: rho1_atom_set, rho_atom_set
206 TYPE(section_vals_type),
POINTER :: input, xc_section, xc_section_aux
208 CALL timeset(routinen, handle)
210 NULLIFY (auxbas_pw_pool, pw_env, v_rspace_new, para_env, rho1_r, &
211 v_xc, rho1_ao, rho_ao, poisson_env, input, rho, dft_control, &
212 logger, rho1_g, v_xc_tau)
215 energy_hartree = 0.0_dp
216 energy_hartree_1c = 0.0_dp
218 cpassert(
ASSOCIATED(p_env%kpp1))
219 cpassert(
ASSOCIATED(p_env%kpp1_env))
220 kpp1_env => p_env%kpp1_env
230 linres_control=linres_control, &
231 dft_control=dft_control)
233 gapw = dft_control%qs_control%gapw
234 gapw_xc = dft_control%qs_control%gapw_xc
235 lr_triplet = linres_control%lr_triplet
238 rho1_xc => p_env%rho1_xc
239 cpassert(
ASSOCIATED(rho1))
241 cpassert(
ASSOCIATED(rho1_xc))
244 CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
245 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
247 nspins =
SIZE(p_env%kpp1)
248 lrigpw = dft_control%qs_control%lrigpw
252 lri_density=lri_density, &
253 atomic_kind_set=atomic_kind_set)
256 IF (.NOT.
ASSOCIATED(kpp1_env%v_ao))
THEN
260 ALLOCATE (kpp1_env%v_ao(ispin)%matrix)
261 CALL dbcsr_copy(kpp1_env%v_ao(ispin)%matrix, matrix_s(1)%matrix, &
262 name=
"kpp1%v_ao-"//adjustl(cp_to_string(ispin)))
266 IF (dft_control%do_admm)
THEN
267 xc_section => admm_env%xc_section_primary
273 cpassert(
ASSOCIATED(pw_env))
274 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
275 poisson_env=poisson_env)
276 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
277 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
279 IF (gapw .OR. gapw_xc) &
283 CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
286 CALL pw_copy(rho1_g(1), rho1_tot_gspace)
288 CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
291 CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
293 IF (.NOT. (nspins == 1 .AND. lr_triplet))
THEN
294 CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
297 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
300 CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
314 IF (deriv2_analytic)
THEN
315 CALL qs_rho_get(rho1a, rho_r=rho1_r, tau_r=tau1_r)
316 CALL qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, lr_triplet, v_xc, v_xc_tau)
317 IF (gapw .OR. gapw_xc)
THEN
318 CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
319 rho1_atom_set => p_env%local_rho_set%rho_atom_set
321 do_tddft=.false., do_triplet=lr_triplet)
324 CALL qs_fxc_fdiff(ks_env, rho0, rho1a, xc_section, 6, lr_triplet, v_xc, v_xc_tau)
325 cpassert((.NOT. gapw) .AND. (.NOT. gapw_xc))
331 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
333 CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
334 IF (
ASSOCIATED(v_xc_tau))
CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
338 IF (dft_control%do_admm)
THEN
340 IF (.NOT.
ASSOCIATED(kpp1_env%deriv_set_admm))
THEN
341 cpassert(.NOT. lr_triplet)
342 xc_section_aux => admm_env%xc_section_aux
345 ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
347 rho_r, auxbas_pw_pool, &
348 xc_section=xc_section_aux)
357 CALL dbcsr_set(kpp1_env%v_ao(ispin)%matrix, 0.0_dp)
363 IF (nspins == 1)
THEN
365 IF (.NOT. (lr_triplet))
THEN
366 CALL pw_scale(v_rspace_new(1), 2.0_dp)
367 IF (
ASSOCIATED(v_xc_tau))
CALL pw_scale(v_xc_tau(1), 2.0_dp)
371 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
372 pmat=rho1_ao(ispin), &
373 hmat=kpp1_env%v_ao(ispin), &
375 calculate_forces=.false., gapw=gapw_xc)
377 IF (
ASSOCIATED(v_xc_tau))
THEN
378 CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
379 pmat=rho1_ao(ispin), &
380 hmat=kpp1_env%v_ao(ispin), &
382 compute_tau=.true., &
383 calculate_forces=.false., gapw=gapw_xc)
387 IF (.NOT. lr_triplet)
THEN
388 CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp, 0.0_dp)
390 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
391 pmat=rho_ao(ispin), &
392 hmat=kpp1_env%v_ao(ispin), &
394 calculate_forces=.false., gapw=gapw)
398 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
399 pmat=rho_ao(ispin), &
400 hmat=kpp1_env%v_ao(ispin), &
402 calculate_forces=.false., gapw=gapw_xc)
404 IF (
ASSOCIATED(v_xc_tau))
THEN
405 CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
406 pmat=rho_ao(ispin), &
407 hmat=kpp1_env%v_ao(ispin), &
409 compute_tau=.true., &
410 calculate_forces=.false., gapw=gapw_xc)
413 CALL pw_copy(v_hartree_rspace, v_rspace_new(ispin))
414 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
415 pmat=rho_ao(ispin), &
416 hmat=kpp1_env%v_ao(ispin), &
418 calculate_forces=.false., gapw=gapw)
423 IF (nspins == 1)
THEN
424 IF (.NOT. (lr_triplet))
THEN
425 CALL pw_scale(v_rspace_new(1), 2.0_dp)
426 IF (
ASSOCIATED(v_xc_tau))
CALL pw_scale(v_xc_tau(1), 2.0_dp)
430 IF (.NOT. lr_triplet)
THEN
431 CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp)
434 CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin), 1.0_dp)
438 IF (
ASSOCIATED(v_xc_tau)) &
439 cpabort(
"metaGGA-functionals not supported with LRI!")
441 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
444 lri_v_int(ikind)%v_int = 0.0_dp
446 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
447 lri_v_int, .false.,
"LRI_AUX")
449 CALL para_env%sum(lri_v_int(ikind)%v_int)
452 k1mat(1)%matrix => kpp1_env%v_ao(ispin)%matrix
453 IF (lri_env%exact_1c_terms)
THEN
454 CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, &
455 rho_ao(ispin)%matrix, qs_env, .false.,
"ORB")
460 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
461 pmat=rho_ao(ispin), &
462 hmat=kpp1_env%v_ao(ispin), &
464 calculate_forces=.false., gapw=gapw)
466 IF (
ASSOCIATED(v_xc_tau))
THEN
467 CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
468 pmat=rho_ao(ispin), &
469 hmat=kpp1_env%v_ao(ispin), &
471 compute_tau=.true., &
472 calculate_forces=.false., gapw=gapw)
478 CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, kpp1_env%v_ao(ispin)%matrix)
482 IF (.NOT. ((nspins == 1 .AND. lr_triplet)))
THEN
484 p_env%hartree_local%ecoul_1c, &
485 p_env%local_rho_set, &
486 para_env, tddft=.true., core_2nd=.true.)
489 calculate_forces=.false., &
490 local_rho_set=p_env%local_rho_set)
494 ns =
SIZE(p_env%kpp1)
495 ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
497 psmat(1:ns, 1:1) => rho_ao(1:ns)
498 CALL update_ks_atom(qs_env, ksmat, psmat, forces=.false., tddft=.true., &
499 rho_atom_external=p_env%local_rho_set%rho_atom_set)
500 ELSEIF (gapw_xc)
THEN
501 ns =
SIZE(p_env%kpp1)
502 ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
504 psmat(1:ns, 1:1) => rho_ao(1:ns)
505 CALL update_ks_atom(qs_env, ksmat, psmat, forces=.false., tddft=.true., &
506 rho_atom_external=p_env%local_rho_set%rho_atom_set)
510 IF (dft_control%qs_control%do_kg .AND. .NOT. (lr_triplet .OR. gapw .OR. gapw_xc))
THEN
519 ks_matrix=p_env%kpp1, &
521 calc_force=.false., &
527 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
528 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
530 CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
532 DEALLOCATE (v_rspace_new)
533 IF (
ASSOCIATED(v_xc_tau))
THEN
535 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
537 DEALLOCATE (v_xc_tau)
540 CALL timestop(handle)
542 END SUBROUTINE apply_op_2_dft
549 SUBROUTINE apply_op_2_xtb(qs_env, p_env)
550 TYPE(qs_environment_type),
POINTER :: qs_env
551 TYPE(qs_p_env_type) :: p_env
553 CHARACTER(len=*),
PARAMETER :: routinen =
'apply_op_2_xtb'
555 INTEGER :: atom_a, handle, iatom, ikind, is, ispin, &
556 na, natom, natorb, nkind, ns, nsgf, &
558 INTEGER,
DIMENSION(25) :: lao
559 INTEGER,
DIMENSION(5) :: occ
560 LOGICAL :: lr_triplet
561 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: mcharge, mcharge1
562 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: aocg, aocg1, charges, charges1
563 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
564 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: rho_ao
565 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_p1, matrix_s
566 TYPE(dft_control_type),
POINTER :: dft_control
567 TYPE(linres_control_type),
POINTER :: linres_control
568 TYPE(mp_para_env_type),
POINTER :: para_env
569 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
570 TYPE(pw_env_type),
POINTER :: pw_env
571 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
572 TYPE(qs_kpp1_env_type),
POINTER :: kpp1_env
573 TYPE(qs_rho_type),
POINTER :: rho, rho1
574 TYPE(xtb_atom_type),
POINTER :: xtb_kind
576 CALL timeset(routinen, handle)
578 cpassert(
ASSOCIATED(p_env%kpp1_env))
579 cpassert(
ASSOCIATED(p_env%kpp1))
580 kpp1_env => p_env%kpp1_env
583 cpassert(
ASSOCIATED(rho1))
589 linres_control=linres_control, &
590 dft_control=dft_control)
594 lr_triplet = linres_control%lr_triplet
595 cpassert(.NOT. lr_triplet)
597 nspins =
SIZE(p_env%kpp1)
600 CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
603 IF (dft_control%qs_control%xtb_control%coulomb_interaction)
THEN
605 CALL get_qs_env(qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
606 natom =
SIZE(particle_set)
609 ALLOCATE (mcharge(natom), charges(natom, 5))
610 ALLOCATE (mcharge1(natom), charges1(natom, 5))
613 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
614 nkind =
SIZE(atomic_kind_set)
616 ALLOCATE (aocg(nsgf, natom))
618 ALLOCATE (aocg1(nsgf, natom))
620 CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
621 CALL ao_charges(matrix_p1, matrix_s, aocg1, para_env)
624 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
627 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
628 charges(atom_a, :) = real(occ(:), kind=
dp)
631 charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
632 charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
636 DEALLOCATE (aocg, aocg1)
638 mcharge(iatom) = sum(charges(iatom, :))
639 mcharge1(iatom) = sum(charges1(iatom, :))
644 DEALLOCATE (charges, mcharge, charges1, mcharge1)
647 CALL timestop(handle)
649 END SUBROUTINE apply_op_2_xtb
659 TYPE(qs_environment_type),
POINTER :: qs_env
660 TYPE(qs_p_env_type) :: p_env
662 CHARACTER(LEN=*),
PARAMETER :: routinen =
'apply_hfx'
664 INTEGER :: handle, ispin, nspins
666 REAL(kind=
dp) :: alpha
667 TYPE(cp_logger_type),
POINTER :: logger
668 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: h1_mat, matrix_s, rho1_ao, work
669 TYPE(dft_control_type),
POINTER :: dft_control
670 TYPE(section_vals_type),
POINTER :: hfx_section, input
672 CALL timeset(routinen, handle)
679 dft_control=dft_control)
680 nspins = dft_control%nspins
687 IF (dft_control%do_admm)
THEN
689 cpabort(
"ADMM: Linear Response needs purification_method=none")
692 cpabort(
"ADMM: Linear Response needs scaling_model=none")
695 cpabort(
"ADMM: Linear Response needs admm_method=basis_projection")
698 rho1_ao => p_env%p1_admm
699 h1_mat => p_env%kpp1_admm
708 ALLOCATE (work(ispin)%matrix)
709 CALL dbcsr_create(work(ispin)%matrix, template=h1_mat(ispin)%matrix)
710 CALL dbcsr_copy(work(ispin)%matrix, h1_mat(ispin)%matrix)
711 CALL dbcsr_set(work(ispin)%matrix, 0.0_dp)
714 CALL hfx_matrix(work, rho1_ao, qs_env, hfx_section)
717 IF (nspins == 2) alpha = 1.0_dp
720 CALL dbcsr_add(h1_mat(ispin)%matrix, work(ispin)%matrix, 1.0_dp, alpha)
727 CALL timestop(handle)
743 SUBROUTINE hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections, external_x_data, ex)
744 TYPE(dbcsr_p_type),
DIMENSION(:),
TARGET :: matrix_ks, rho_ao
745 TYPE(qs_environment_type),
POINTER :: qs_env
746 TYPE(section_vals_type),
POINTER :: hfx_sections
747 TYPE(hfx_type),
DIMENSION(:, :),
OPTIONAL,
TARGET :: external_x_data
748 REAL(kind=
dp),
OPTIONAL :: ex
750 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_matrix'
752 INTEGER :: handle, irep, ispin, mspin, n_rep_hf, &
754 LOGICAL :: distribute_fock_matrix, &
755 hfx_treat_lsd_in_core, &
757 REAL(kind=
dp) :: eh1, ehfx
758 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks_kp, rho_ao_kp
759 TYPE(dft_control_type),
POINTER :: dft_control
760 TYPE(hfx_type),
DIMENSION(:, :),
POINTER :: x_data
761 TYPE(mp_para_env_type),
POINTER :: para_env
763 CALL timeset(routinen, handle)
765 NULLIFY (dft_control, para_env, matrix_ks_kp, rho_ao_kp, x_data)
768 dft_control=dft_control, &
770 s_mstruct_changed=s_mstruct_changed, &
773 IF (
PRESENT(external_x_data)) x_data => external_x_data
775 cpassert(dft_control%nimages == 1)
776 nspins = dft_control%nspins
783 distribute_fock_matrix = .true.
786 IF (hfx_treat_lsd_in_core) mspin = nspins
788 matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
789 rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
791 DO irep = 1, n_rep_hf
794 IF (x_data(irep, 1)%do_hfx_ri)
THEN
795 CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, &
796 rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
797 nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
803 s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
811 IF (
PRESENT(ex)) ex = ehfx
813 CALL timestop(handle)
823 TYPE(qs_environment_type),
POINTER :: qs_env
824 TYPE(qs_p_env_type) :: p_env
826 CHARACTER(len=*),
PARAMETER :: routinen =
'apply_xc_admm'
828 CHARACTER(LEN=default_string_length) :: basis_type
829 INTEGER :: handle, ispin, ns, nspins
830 INTEGER,
DIMENSION(2, 3) :: bo
832 REAL(kind=
dp) :: alpha
833 TYPE(admm_type),
POINTER :: admm_env
834 TYPE(dbcsr_p_type) :: xcmat
835 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
836 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, psmat
837 TYPE(dft_control_type),
POINTER :: dft_control
838 TYPE(linres_control_type),
POINTER :: linres_control
839 TYPE(mp_para_env_type),
POINTER :: para_env
840 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
841 POINTER :: sab_aux_fit
842 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho1_aux_g
843 TYPE(pw_env_type),
POINTER :: pw_env
844 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
845 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho1_aux_r, tau_pw, v_xc, v_xc_tau
846 TYPE(rho_atom_type),
DIMENSION(:),
POINTER :: rho1_atom_set, rho_atom_set
847 TYPE(section_vals_type),
POINTER :: xc_fun_section, xc_section
848 TYPE(task_list_type),
POINTER :: task_list
849 TYPE(xc_rho_cflags_type) :: needs
850 TYPE(xc_rho_set_type) :: rho1_set
852 CALL timeset(routinen, handle)
854 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
856 IF (dft_control%do_admm)
THEN
860 CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
861 cpassert(.NOT. dft_control%qs_control%lrigpw)
862 cpassert(.NOT. linres_control%lr_triplet)
864 nspins = dft_control%nspins
868 cpassert(
ASSOCIATED(pw_env))
869 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
873 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s)
874 ALLOCATE (xcmat%matrix)
875 CALL dbcsr_create(xcmat%matrix, template=matrix_s(1)%matrix)
878 gapw = admm_env%do_gapw
880 CALL qs_rho_get(p_env%rho1_admm, rho_r=rho1_aux_r, rho_g=rho1_aux_g)
881 xc_section => admm_env%xc_section_aux
882 bo = rho1_aux_r(1)%pw_grid%bounds_local
897 CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, p_env%kpp1_env%deriv_set_admm, p_env%kpp1_env%rho_set_admm, &
898 rho1_aux_r, rho1_aux_g, tau_pw, auxbas_pw_pool, gapw=.false., &
899 xc_section=xc_section)
900 IF (
ASSOCIATED(v_xc_tau))
THEN
901 cpabort(
"Meta-GGA ADMM functionals not yet supported!")
905 basis_type =
"AUX_FIT"
907 CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
908 IF (admm_env%do_gapw)
THEN
910 do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
911 rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
912 rho1_atom_set => p_env%local_rho_set_admm%rho_atom_set
914 kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
915 basis_type =
"AUX_FIT_SOFT"
916 task_list => admm_env%admm_gapw_env%task_list
920 IF (nspins == 1) alpha = 2.0_dp
923 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
924 CALL dbcsr_copy(xcmat%matrix, matrix_s(1)%matrix)
925 CALL dbcsr_set(xcmat%matrix, 0.0_dp)
926 CALL integrate_v_rspace(v_rspace=v_xc(ispin), hmat=xcmat, qs_env=qs_env, &
927 calculate_forces=.false., basis_type=basis_type, &
928 task_list_external=task_list)
929 CALL dbcsr_add(p_env%kpp1_admm(ispin)%matrix, xcmat%matrix, 1.0_dp, alpha)
932 IF (admm_env%do_gapw)
THEN
934 ns =
SIZE(p_env%kpp1_admm)
935 ksmat(1:ns, 1:1) => p_env%kpp1_admm(1:ns)
936 psmat(1:ns, 1:1) => p_env%p1_admm(1:ns)
937 CALL update_ks_atom(qs_env, ksmat, psmat, forces=.false., tddft=.true., &
938 rho_atom_external=p_env%local_rho_set_admm%rho_atom_set, &
939 kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
940 oce_external=admm_env%admm_gapw_env%oce, &
941 sab_external=sab_aux_fit)
945 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
948 CALL dbcsr_deallocate_matrix(xcmat%matrix)
953 CALL timestop(handle)
Types and set/get functions for auxiliary density matrix methods.
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
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
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
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
subroutine, public vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, core_2nd)
Calculates one center GAPW Hartree energies and matrix elements Hartree potentials are input Takes po...
Routines to calculate HFX energy and potential.
subroutine, public integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, geometry_did_change, irep, distribute_fock_matrix, ispin)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
subroutine, public hfx_ri_update_ks(qs_env, ri_data, ks_matrix, ehfx, mos, rho_ao, geometry_did_change, nspins, hf_fraction)
...
Types and set/get functions for HFX.
Routines for a Kim-Gordon-like partitioning into molecular subunits.
subroutine, public kg_ekin_subset(qs_env, ks_matrix, ekin_mol, calc_force, do_kernel, pmat_ext)
Calculates the subsystem Hohenberg-Kohn kinetic energy and the forces.
Types needed for a Kim-Gordon-like partitioning into molecular subunits.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
routines that build the Kohn-Sham matrix for the LRIGPW and xc parts
subroutine, public calculate_lri_ks_matrix(lri_env, lri_v_int, h_matrix, atomic_kind_set, cell_to_index)
update of LRIGPW KS matrix
Interface to the message passing library MPI.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Define the data structure for the particle information.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
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.
https://en.wikipedia.org/wiki/Finite_difference_coefficient
subroutine, public qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau)
...
subroutine, public qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, fxc_rho, fxc_tau)
...
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
...
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
basis types for the calculation of the perturbation of density theory.
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
subroutine, public apply_xc_admm(qs_env, p_env)
...
subroutine, public hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections, external_x_data, ex)
Add the hfx contributions to the Hamiltonian.
subroutine, public apply_hfx(qs_env, p_env)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public apply_op_2(qs_env, p_env, c0, Av)
...
Type definitiona for linear response calculations.
Define the neighbor list data types and the corresponding functionality.
Utility functions for the perturbation calculations.
subroutine, public p_env_finish_kpp1(qs_env, p_env)
...
basis types for the calculation of the perturbation of density theory.
subroutine, public integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, local_rho_set_2nd, atener, kforce)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
subroutine, public calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, do_tddft, do_tddfpt2, do_triplet, kind_set_external)
...
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Exchange and Correlation functional calculations.
subroutine, public xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, pw_pool, xc_section, gapw, vxg, lsd_singlets, do_excitations, do_triplet, do_tddft, compute_virial, virial_xc)
Caller routine to calculate the second order potential in the direction of rho1_r.
subroutine, public xc_prep_2nd_deriv(deriv_set, rho_set, rho_r, pw_pool, xc_section, tau_r)
Prepare objects for the calculation of the 2nd derivatives of the density functional....
Calculation of Coulomb Hessian contributions in xTB.
subroutine, public xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge)
...
Definition of the xTB parameter types.
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, occupation, electronegativity, chmax)
...