108 integrate_v_rspace_diagonal,&
109 integrate_v_rspace_one_center
123#include "./base/base_uses.f90"
129 LOGICAL,
PARAMETER :: debug_this_module = .true.
130 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_ks_utils'
148 SUBROUTINE low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
149 calculate_forces, auxbas_pw_pool)
154 LOGICAL,
INTENT(IN) :: do_hfx, just_energy, calculate_forces
157 CHARACTER(*),
PARAMETER :: routinen =
'low_spin_roks'
159 INTEGER :: handle, irep, ispin, iterm, k, k_alpha, &
160 k_beta, n_rep, nelectron, nspin, nterms
161 INTEGER,
DIMENSION(:),
POINTER :: ivec
162 INTEGER,
DIMENSION(:, :, :),
POINTER :: occupations
163 LOGICAL :: compute_virial, in_range, &
165 REAL(kind=
dp) :: ehfx, exc
166 REAL(kind=
dp),
DIMENSION(3, 3) :: virial_xc_tmp
167 REAL(kind=
dp),
DIMENSION(:),
POINTER :: energy_scaling, rvec, scaling
169 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_h, matrix_hfx, matrix_p, mdummy, &
171 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p2
172 TYPE(
dbcsr_type),
POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, &
174 TYPE(
hfx_type),
DIMENSION(:, :),
POINTER :: x_data
175 TYPE(
mo_set_type),
DIMENSION(:),
POINTER :: mo_array
181 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, tau, vxc, vxc_tau
185 low_spin_roks_section, xc_section
188 IF (.NOT. dft_control%low_spin_roks)
RETURN
190 CALL timeset(routinen, handle)
192 NULLIFY (ks_env, rho_ao)
195 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc)
THEN
196 CALL cp_abort(__location__,
"GAPW/GAPW_XC are not compatible with low spin ROKS method.")
198 IF (dft_control%do_admm)
THEN
199 CALL cp_abort(__location__,
"ADMM not compatible with low spin ROKS method.")
201 IF (dft_control%do_admm)
THEN
203 CALL cp_abort(__location__,
"ADMM with XC correction functional "// &
204 "not compatible with low spin ROKS method.")
207 IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
208 dft_control%qs_control%xtb)
THEN
209 CALL cp_abort(__location__,
"SE/xTB/DFTB are not compatible with low spin ROKS method.")
214 mo_derivs=mo_derivs, &
224 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
230 cpassert(
SIZE(mo_array, 1) == 2)
233 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
234 cpassert(uniform_occupation)
235 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
236 cpassert(uniform_occupation)
237 IF (do_hfx .AND. calculate_forces .AND. compute_virial)
THEN
238 CALL cp_abort(__location__,
"ROKS virial with HFX not available.")
241 NULLIFY (dbcsr_deriv)
243 CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix)
247 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
249 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff)
257 ALLOCATE (energy_scaling(nterms))
258 energy_scaling = rvec
261 cpassert(n_rep == nterms)
262 CALL section_vals_val_get(low_spin_roks_section,
"SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec)
263 nelectron =
SIZE(ivec)
264 cpassert(nelectron == k_alpha - k_beta)
265 ALLOCATE (occupations(2, nelectron, nterms))
268 CALL section_vals_val_get(low_spin_roks_section,
"SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec)
269 cpassert(nelectron ==
SIZE(ivec))
270 in_range = all(ivec >= 1) .AND. all(ivec <= 2)
273 occupations(ivec(k), k, iterm) = 1
283 ALLOCATE (matrix_p(ispin)%matrix)
284 CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, &
285 name=
"density matrix low spin roks")
286 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
292 ALLOCATE (matrix_h(ispin)%matrix)
293 CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, &
294 name=
"KS matrix low spin roks")
295 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
302 ALLOCATE (matrix_hfx(ispin)%matrix)
303 CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
304 name=
"HFX matrix low spin roks")
310 NULLIFY (tau, vxc_tau, vxc)
311 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
313 ALLOCATE (rho_r(nspin))
314 ALLOCATE (rho_g(nspin))
316 CALL auxbas_pw_pool%create_pw(rho_r(ispin))
317 CALL auxbas_pw_pool%create_pw(rho_g(ispin))
319 CALL auxbas_pw_pool%create_pw(work_v_rspace)
323 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
324 NULLIFY (fm_scaled, fm_deriv)
330 ALLOCATE (scaling(k_alpha))
337 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
339 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
343 0.0_dp, matrix_p(ispin)%matrix, retain_sparsity=.true.)
346 rho=rho_r(ispin), rho_gspace=rho_g(ispin), &
351 IF (just_energy)
THEN
352 exc =
xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
355 cpassert(.NOT. compute_virial)
357 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
358 pw_pool=xc_pw_pool, compute_virial=.false., virial_xc=virial_xc_tmp)
361 energy%exc = energy%exc + energy_scaling(iterm)*exc
366 CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
370 recalc_integrals=.false., update_energy=.true.)
371 energy%ex = ehfx + energy_scaling(iterm)*energy%ex
375 IF (.NOT. just_energy)
THEN
378 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
380 CALL pw_axpy(vxc(ispin), work_v_rspace, energy_scaling(iterm)*vxc(ispin)%pw_grid%dvol, 0.0_dp)
381 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
382 qs_env=qs_env, calculate_forces=calculate_forces)
383 CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
390 CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
391 1.0_dp, energy_scaling(iterm))
393 IF (calculate_forces)
THEN
394 CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
395 IF (x_data(1, 1)%n_rep_hf /= 1)
THEN
396 CALL cp_abort(__location__,
"Multiple HFX section forces not compatible "// &
397 "with low spin ROKS method.")
399 IF (x_data(1, 1)%do_hfx_ri)
THEN
400 CALL cp_abort(__location__,
"HFX_RI forces not compatible with low spin ROKS method.")
404 matrix_p2(1:nspin, 1:1) => matrix_p(1:nspin)
406 irep, compute_virial, &
407 adiabatic_rescale_factor=energy_scaling(iterm))
415 CALL dbcsr_multiply(
'n',
'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
416 0.0_dp, dbcsr_deriv, last_column=k_alpha)
419 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
421 CALL dbcsr_add(mo_derivs(1)%matrix, dbcsr_deriv, 1.0_dp, 1.0_dp)
430 CALL auxbas_pw_pool%give_back_pw(rho_r(ispin))
431 CALL auxbas_pw_pool%give_back_pw(rho_g(ispin))
433 DEALLOCATE (rho_r, rho_g)
440 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
445 DEALLOCATE (occupations)
446 DEALLOCATE (energy_scaling)
451 CALL timestop(handle)
465 calculate_forces, auxbas_pw_pool)
471 LOGICAL,
INTENT(IN) :: just_energy, calculate_forces
474 CHARACTER(*),
PARAMETER :: routinen =
'sic_explicit_orbitals'
476 INTEGER :: handle, i, iorb, k_alpha, k_beta, norb
477 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: sic_orbital_list
478 LOGICAL :: compute_virial, uniform_occupation
479 REAL(kind=
dp) :: ener, exc
480 REAL(kind=
dp),
DIMENSION(3, 3) :: virial_xc_tmp
484 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: mo_derivs_local
487 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: mo_derivs, rho_ao, tmp_dbcsr
488 TYPE(
dbcsr_type),
POINTER :: orb_density_matrix, orb_h
489 TYPE(
mo_set_type),
DIMENSION(:),
POINTER :: mo_array
496 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, tau, vxc, vxc_tau
503 IF (dft_control%sic_method_id .NE.
sic_eo)
RETURN
505 CALL timeset(routinen, handle)
507 NULLIFY (tau, vxc_tau, mo_derivs, ks_env, rho_ao)
512 mo_derivs=mo_derivs, &
522 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
525 DO i = 1,
SIZE(mo_array)
526 IF (mo_array(i)%use_mo_coeff_b)
THEN
528 mo_array(i)%mo_coeff)
532 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
535 cpassert(
SIZE(mo_array, 1) == 2)
537 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
538 cpassert(uniform_occupation)
539 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff, uniform_occupation=uniform_occupation)
540 cpassert(uniform_occupation)
544 DO i = 1,
SIZE(mo_derivs, 1)
546 NULLIFY (tmp_dbcsr(i)%matrix)
548 CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix)
549 CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp)
552 k_alpha = 0; k_beta = 0
553 SELECT CASE (dft_control%sic_list_id)
556 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
559 IF (
SIZE(mo_array, 1) > 1)
THEN
560 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
564 norb = k_alpha + k_beta
565 ALLOCATE (sic_orbital_list(3, norb))
570 sic_orbital_list(1, iorb) = 1
571 sic_orbital_list(2, iorb) = i
572 sic_orbital_list(3, iorb) = 1
576 sic_orbital_list(1, iorb) = 2
577 sic_orbital_list(2, iorb) = i
578 IF (
SIZE(mo_derivs, 1) == 1)
THEN
579 sic_orbital_list(3, iorb) = 1
581 sic_orbital_list(3, iorb) = 2
587 cpassert(
SIZE(mo_array, 1) == 2)
589 cpassert(
SIZE(mo_derivs, 1) == 1)
590 cpassert(dft_control%restricted)
592 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
595 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
598 norb = k_alpha - k_beta
599 ALLOCATE (sic_orbital_list(3, norb))
602 DO i = k_beta + 1, k_alpha
604 sic_orbital_list(1, iorb) = 1
605 sic_orbital_list(2, iorb) = i
607 sic_orbital_list(3, iorb) = 1
615 CALL auxbas_pw_pool%create_pw(orb_rho_r)
616 CALL auxbas_pw_pool%create_pw(tmp_r)
617 CALL auxbas_pw_pool%create_pw(orb_rho_g)
618 CALL auxbas_pw_pool%create_pw(tmp_g)
619 CALL auxbas_pw_pool%create_pw(work_v_gspace)
620 CALL auxbas_pw_pool%create_pw(work_v_rspace)
622 ALLOCATE (orb_density_matrix)
623 CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, &
624 name=
"orb_density_matrix")
625 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
626 orb_density_matrix_p%matrix => orb_density_matrix
630 name=
"orb_density_matrix")
632 orb_h_p%matrix => orb_h
634 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
637 template_fmstruct=mo_coeff%matrix_struct)
638 CALL cp_fm_create(matrix_v, fm_struct_tmp, name=
"matrix_v")
639 CALL cp_fm_create(matrix_hv, fm_struct_tmp, name=
"matrix_hv")
642 ALLOCATE (mo_derivs_local(
SIZE(mo_array, 1)))
643 DO i = 1,
SIZE(mo_array, 1)
644 CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff)
645 CALL cp_fm_create(mo_derivs_local(i), mo_coeff%matrix_struct)
662 CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb)), mo_coeff=mo_coeff)
663 CALL cp_fm_to_fm(mo_coeff, matrix_v, 1, sic_orbital_list(2, iorb), 1)
666 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
671 rho=orb_rho_r, rho_gspace=orb_rho_g, &
679 energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener
680 IF (.NOT. just_energy)
THEN
682 CALL pw_scale(work_v_rspace, -dft_control%sic_scaling_a*work_v_rspace%pw_grid%dvol)
686 IF (just_energy)
THEN
687 exc =
xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
690 cpassert(.NOT. compute_virial)
692 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
693 pw_pool=xc_pw_pool, compute_virial=compute_virial, virial_xc=virial_xc_tmp)
695 CALL pw_axpy(vxc(1), work_v_rspace, -dft_control%sic_scaling_b*vxc(1)%pw_grid%dvol)
697 energy%exc = energy%exc - dft_control%sic_scaling_b*exc
699 IF (.NOT. just_energy)
THEN
701 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=orb_density_matrix_p, hmat=orb_h_p, &
702 qs_env=qs_env, calculate_forces=calculate_forces)
707 CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb)), 0.0_dp)
708 CALL cp_fm_to_fm(matrix_hv, mo_derivs_local(sic_orbital_list(3, iorb)), 1, 1, sic_orbital_list(2, iorb))
710 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix)
711 CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, &
712 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp)
715 CALL xc_pw_pool%give_back_pw(vxc(1))
716 CALL xc_pw_pool%give_back_pw(vxc(2))
723 CALL auxbas_pw_pool%give_back_pw(orb_rho_r)
724 CALL auxbas_pw_pool%give_back_pw(tmp_r)
725 CALL auxbas_pw_pool%give_back_pw(orb_rho_g)
726 CALL auxbas_pw_pool%give_back_pw(tmp_g)
727 CALL auxbas_pw_pool%give_back_pw(work_v_gspace)
728 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
740 CALL timestop(handle)
757 qs_env, dft_control, rho, poisson_env, just_energy, &
758 calculate_forces, auxbas_pw_pool)
766 LOGICAL,
INTENT(IN) :: just_energy, calculate_forces
769 INTEGER :: i, nelec, nelec_a, nelec_b, nforce
770 REAL(kind=
dp) :: ener, full_scaling, scaling
771 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: store_forces
772 TYPE(
mo_set_type),
DIMENSION(:),
POINTER :: mo_array
777 NULLIFY (mo_array, rho_g)
779 IF (dft_control%sic_method_id ==
sic_none)
RETURN
780 IF (dft_control%sic_method_id ==
sic_eo)
RETURN
782 IF (dft_control%qs_control%gapw) &
783 cpabort(
"sic and GAPW not yet compatible")
786 cpassert(dft_control%nspins == 2)
788 CALL auxbas_pw_pool%create_pw(work_rho)
789 CALL auxbas_pw_pool%create_pw(work_v)
794 SELECT CASE (dft_control%sic_method_id)
796 CALL pw_copy(rho_g(1), work_rho)
797 CALL pw_axpy(rho_g(2), work_rho, alpha=-1._dp)
802 CALL get_mo_set(mo_set=mo_array(1), nelectron=nelec_a)
803 CALL get_mo_set(mo_set=mo_array(2), nelectron=nelec_b)
804 nelec = nelec_a + nelec_b
805 CALL pw_copy(rho_g(1), work_rho)
806 CALL pw_axpy(rho_g(2), work_rho)
807 scaling = 1.0_dp/real(nelec, kind=
dp)
811 cpabort(
"Unknown sic method id")
816 IF (calculate_forces)
THEN
819 DO i = 1,
SIZE(force)
820 nforce = nforce +
SIZE(force(i)%ch_pulay, 2)
822 ALLOCATE (store_forces(3, nforce))
824 DO i = 1,
SIZE(force)
825 store_forces(1:3, nforce + 1:nforce +
SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :)
826 force(i)%ch_pulay(:, :) = 0.0_dp
827 nforce = nforce +
SIZE(force(i)%ch_pulay, 2)
834 v_hartree_gspace=work_v, &
835 calculate_forces=calculate_forces, &
836 itype_of_density=
"SPIN")
838 SELECT CASE (dft_control%sic_method_id)
840 full_scaling = -dft_control%sic_scaling_a
842 full_scaling = -dft_control%sic_scaling_a*nelec
844 cpabort(
"Unknown sic method id")
846 energy%hartree = energy%hartree + full_scaling*ener
849 IF (calculate_forces)
THEN
851 DO i = 1,
SIZE(force)
852 force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + &
853 store_forces(1:3, nforce + 1:nforce +
SIZE(force(i)%ch_pulay, 2))
854 nforce = nforce +
SIZE(force(i)%ch_pulay, 2)
858 IF (.NOT. just_energy)
THEN
859 ALLOCATE (v_sic_rspace)
860 CALL auxbas_pw_pool%create_pw(v_sic_rspace)
864 dft_control%sic_scaling_a*v_sic_rspace%pw_grid%dvol)
867 CALL auxbas_pw_pool%give_back_pw(work_rho)
868 CALL auxbas_pw_pool%give_back_pw(work_v)
881 INTEGER :: img, ispin, n_electrons, output_unit
882 REAL(
dp) :: tot1_h, tot1_s, tot_rho_r, trace, &
884 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tot_rho_r_arr
887 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_s, rho_ao
890 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
893 NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
894 dft_control, tot_rho_r_arr, rho_ao)
899 qs_kind_set=qs_kind_set, &
900 cell=cell, qs_charges=qs_charges, &
902 matrix_s_kp=matrix_s, &
903 dft_control=dft_control)
911 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
912 n_electrons = n_electrons - dft_control%charge
917 DO ispin = 1, dft_control%nspins
918 DO img = 1, dft_control%nimages
919 CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
920 trace = trace + trace_tmp
925 IF (output_unit > 0)
THEN
926 WRITE (unit=output_unit, fmt=
"(/,T3,A,T41,F20.10)")
"Trace(PS):", trace
927 WRITE (unit=output_unit, fmt=
"((T3,A,T41,2F20.10))") &
928 "Electronic density on regular grids: ", &
931 REAL(n_electrons,
dp), &
932 "Core density on regular grids:", &
933 qs_charges%total_rho_core_rspace, &
934 qs_charges%total_rho_core_rspace + &
935 qs_charges%total_rho1_hard_nuc - &
936 REAL(n_electrons + dft_control%charge,
dp)
938 IF (dft_control%qs_control%gapw)
THEN
939 tot1_h = qs_charges%total_rho1_hard(1)
940 tot1_s = qs_charges%total_rho1_soft(1)
941 DO ispin = 2, dft_control%nspins
942 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
943 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
945 IF (output_unit > 0)
THEN
946 WRITE (unit=output_unit, fmt=
"((T3,A,T41,2F20.10))") &
947 "Hard and soft densities (Lebedev):", &
949 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
950 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
951 tot_rho_r + tot1_h - tot1_s, &
952 "Total charge density (r-space): ", &
953 tot_rho_r + tot1_h - tot1_s &
954 + qs_charges%total_rho_core_rspace &
955 + qs_charges%total_rho1_hard_nuc
956 IF (qs_charges%total_rho1_hard_nuc /= 0.0_dp)
THEN
957 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
958 "Total CNEO nuc. char. den. (Lebedev): ", &
959 qs_charges%total_rho1_hard_nuc, &
960 "Total CNEO soft char. den. (Lebedev): ", &
961 qs_charges%total_rho1_soft_nuc_lebedev, &
962 "Total CNEO soft char. den. (r-space): ", &
963 qs_charges%total_rho1_soft_nuc_rspace, &
964 "Total soft Rho_e+n+0 (g-space):", &
965 qs_charges%total_rho_gspace
967 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
968 "Total Rho_soft + Rho0_soft (g-space):", &
969 qs_charges%total_rho_gspace
972 qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
973 qs_charges%total_rho_core_rspace + &
974 qs_charges%total_rho1_hard_nuc
976 ELSE IF (dft_control%qs_control%gapw_xc)
THEN
977 tot1_h = qs_charges%total_rho1_hard(1)
978 tot1_s = qs_charges%total_rho1_soft(1)
979 DO ispin = 2, dft_control%nspins
980 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
981 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
983 IF (output_unit > 0)
THEN
984 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T41,2F20.10))") &
985 "Hard and soft densities (Lebedev):", &
987 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
988 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
991 qs_charges%background = tot_rho_r + &
992 qs_charges%total_rho_core_rspace
994 IF (output_unit > 0)
THEN
995 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
996 "Total charge density on r-space grids: ", &
998 qs_charges%total_rho_core_rspace, &
999 "Total charge density g-space grids: ", &
1000 qs_charges%total_rho_gspace
1002 qs_charges%background = tot_rho_r + &
1003 qs_charges%total_rho_core_rspace
1005 IF (output_unit > 0)
WRITE (unit=output_unit, fmt=
"()")
1006 qs_charges%background = qs_charges%background/cell%deth
1009 "PRINT%TOTAL_DENSITIES")
1031 REAL(kind=
dp),
INTENT(IN) :: mulliken_order_p
1033 INTEGER :: bc, n, output_unit, psolver
1034 REAL(kind=
dp) :: ddapc_order_p, implicit_ps_ehartree, &
1042 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1043 psolver = pw_env%poisson_env%parameters%solver
1046 extension=
".scfLog")
1047 IF (output_unit > 0)
THEN
1048 IF (dft_control%do_admm)
THEN
1049 WRITE (unit=output_unit, fmt=
"((T3,A,T60,F20.10))") &
1050 "Wfn fit exchange-correlation energy: ", energy%exc_aux_fit
1051 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc)
THEN
1052 WRITE (unit=output_unit, fmt=
"((T3,A,T60,F20.10))") &
1053 "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit
1056 IF (dft_control%do_admm)
THEN
1058 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1059 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1062 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1063 "Core Hamiltonian energy: ", energy%core, &
1064 "Hartree energy: ", implicit_ps_ehartree, &
1065 "Electric enthalpy: ", energy%hartree, &
1066 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1068 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1069 "Core Hamiltonian energy: ", energy%core, &
1070 "Hartree energy: ", energy%hartree, &
1071 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1074 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1075 "Core Hamiltonian energy: ", energy%core, &
1076 "Hartree energy: ", energy%hartree, &
1077 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1081 IF (dft_control%apply_external_density)
THEN
1082 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1083 "DOING ZMP CALCULATION FROM EXTERNAL DENSITY "
1084 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1085 "Core Hamiltonian energy: ", energy%core, &
1086 "Hartree energy: ", energy%hartree
1087 ELSE IF (dft_control%apply_external_vxc)
THEN
1088 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1089 "DOING ZMP READING EXTERNAL VXC "
1090 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1091 "Core Hamiltonian energy: ", energy%core, &
1092 "Hartree energy: ", energy%hartree
1095 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1096 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1099 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1100 "Core Hamiltonian energy: ", energy%core, &
1101 "Hartree energy: ", implicit_ps_ehartree, &
1102 "Electric enthalpy: ", energy%hartree, &
1103 "Exchange-correlation energy: ", energy%exc
1105 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1106 "Core Hamiltonian energy: ", energy%core, &
1107 "Hartree energy: ", energy%hartree, &
1108 "Exchange-correlation energy: ", energy%exc
1111 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1112 "Core Hamiltonian energy: ", energy%core, &
1113 "Hartree energy: ", energy%hartree, &
1114 "Exchange-correlation energy: ", energy%exc
1119 IF (dft_control%apply_external_density)
THEN
1120 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1121 "Integral of the (density * v_xc): ", energy%exc
1124 IF (energy%e_hartree /= 0.0_dp) &
1125 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1126 "Coulomb (electron-electron) energy: ", energy%e_hartree
1127 IF (energy%dispersion /= 0.0_dp) &
1128 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1129 "Dispersion energy: ", energy%dispersion
1130 IF (energy%efield /= 0.0_dp) &
1131 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1132 "Electric field interaction energy: ", energy%efield
1133 IF (energy%gcp /= 0.0_dp) &
1134 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1135 "gCP energy: ", energy%gcp
1136 IF (dft_control%qs_control%gapw)
THEN
1137 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1138 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit, &
1139 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
1141 IF (dft_control%qs_control%gapw_xc)
THEN
1142 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T61,F20.10))") &
1143 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit
1145 IF (dft_control%dft_plus_u)
THEN
1146 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1147 "DFT+U energy:", energy%dft_plus_u
1149 IF (qs_env%qmmm)
THEN
1150 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1151 "QM/MM Electrostatic energy: ", energy%qmmm_el
1152 IF (qs_env%qmmm_env_qm%image_charge)
THEN
1153 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1154 "QM/MM image charge energy: ", energy%image_charge
1157 IF (dft_control%qs_control%mulliken_restraint)
THEN
1158 WRITE (unit=output_unit, fmt=
"(T3,A,T41,2F20.10)") &
1159 "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken
1161 IF (dft_control%qs_control%ddapc_restraint)
THEN
1162 DO n = 1,
SIZE(dft_control%qs_control%ddapc_restraint_control)
1164 dft_control%qs_control%ddapc_restraint_control(n)%ddapc_order_p
1165 WRITE (unit=output_unit, fmt=
"(T3,A,T41,2F20.10)") &
1166 "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n)
1169 IF (dft_control%qs_control%s2_restraint)
THEN
1170 s2_order_p = dft_control%qs_control%s2_restraint_control%s2_order_p
1171 WRITE (unit=output_unit, fmt=
"(T3,A,T41,2F20.10)") &
1172 "S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint
1174 IF (energy%core_cneo /= 0.0_dp)
THEN
1175 WRITE (unit=output_unit, fmt=
"(T3,A,T61,F20.10)") &
1176 "CNEO| quantum nuclear core energy: ", energy%core_cneo
1181 "DFT%SCF%PRINT%DETAILED_ENERGY")
1199 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_vxc
1201 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_matrix_vxc'
1203 INTEGER :: handle, ispin
1208 CALL timeset(routinen, handle)
1211 IF (
ASSOCIATED(matrix_vxc))
THEN
1215 ALLOCATE (matrix_vxc(
SIZE(matrix_ks)))
1216 DO ispin = 1,
SIZE(matrix_ks)
1217 NULLIFY (matrix_vxc(ispin)%matrix)
1219 CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
1221 CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
1225 CALL get_qs_env(qs_env, dft_control=dft_control)
1226 gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1227 DO ispin = 1,
SIZE(matrix_ks)
1228 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1229 hmat=matrix_vxc(ispin), &
1231 calculate_forces=.false., &
1234 CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
1237 CALL timestop(handle)
1265 vppl_rspace, v_rspace_new, &
1266 v_rspace_new_aux_fit, v_tau_rspace, &
1267 v_tau_rspace_aux_fit, &
1268 v_sic_rspace, v_spin_ddapc_rest_r, &
1269 v_sccs_rspace, v_rspace_embed, cdft_control, &
1273 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
1277 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
POINTER :: v_rspace_new, v_rspace_new_aux_fit, &
1278 v_tau_rspace, v_tau_rspace_aux_fit
1279 TYPE(
pw_r3d_rs_type),
POINTER :: v_sic_rspace, v_spin_ddapc_rest_r, &
1283 LOGICAL,
INTENT(in) :: calculate_forces
1285 CHARACTER(LEN=*),
PARAMETER :: routinen =
'sum_up_and_integrate'
1287 CHARACTER(LEN=default_string_length) :: basis_type
1288 INTEGER :: handle, igroup, ikind, img, ispin, &
1290 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1291 LOGICAL :: do_ppl, gapw, gapw_xc, lrigpw, rigpw
1292 REAL(kind=
dp) :: csign, dvol, fadm
1295 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: ksmat, rho_ao, rho_ao_nokp, smat
1296 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks_aux_fit, &
1297 matrix_ks_aux_fit_dft, rho_ao_aux, &
1313 CALL timeset(routinen, handle)
1315 NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
1316 v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
1317 ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
1318 rho_ao_nokp, ks_env, admm_env, task_list)
1321 dft_control=dft_control, &
1323 v_hartree_rspace=v_rspace, &
1327 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
1328 gapw = dft_control%qs_control%gapw
1329 gapw_xc = dft_control%qs_control%gapw_xc
1330 do_ppl = dft_control%qs_control%do_ppl_method ==
do_ppl_grid
1332 rigpw = dft_control%qs_control%rigpw
1333 lrigpw = dft_control%qs_control%lrigpw
1334 IF (lrigpw .OR. rigpw)
THEN
1337 lri_density=lri_density, &
1338 atomic_kind_set=atomic_kind_set)
1341 nspins = dft_control%nspins
1344 IF (
ASSOCIATED(v_rspace_new))
THEN
1345 DO ispin = 1, nspins
1348 cpassert(dft_control%sic_method_id ==
sic_none)
1350 CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
1353 rho_ao => rho_ao_kp(ispin, :)
1354 ksmat => ks_matrix(ispin, :)
1355 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1356 pmat_kp=rho_ao, hmat_kp=ksmat, &
1358 calculate_forces=calculate_forces, &
1362 CALL pw_copy(v_rspace, v_rspace_new(ispin))
1365 CALL pw_axpy(v_rspace, v_rspace_new(ispin), 1.0_dp, v_rspace_new(ispin)%pw_grid%dvol)
1367 IF (dft_control%qs_control%ddapc_explicit_potential)
THEN
1368 IF (dft_control%qs_control%ddapc_restraint_is_spin)
THEN
1369 IF (ispin == 1)
THEN
1370 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1372 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), -1.0_dp)
1375 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1379 IF (dft_control%qs_control%cdft)
THEN
1380 DO igroup = 1,
SIZE(cdft_control%group)
1381 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1385 IF (ispin == 1)
THEN
1392 IF (ispin == 2) cycle
1395 IF (ispin == 1) cycle
1397 cpabort(
"Unknown constraint type.")
1399 CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_new(ispin), &
1400 csign*cdft_control%strength(igroup))
1407 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1408 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_new(ispin), dvol)
1411 IF (dft_control%do_sccs)
THEN
1412 CALL pw_axpy(v_sccs_rspace, v_rspace_new(ispin))
1415 IF (dft_control%apply_external_potential)
THEN
1417 v_qmmm=vee, scale=-1.0_dp)
1420 cpassert(.NOT. gapw)
1421 CALL pw_axpy(vppl_rspace, v_rspace_new(ispin), vppl_rspace%pw_grid%dvol)
1424 SELECT CASE (dft_control%sic_method_id)
1428 IF (ispin == 1)
THEN
1429 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1431 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), 1.0_dp)
1434 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1439 IF (dft_control%apply_embed_pot)
THEN
1440 CALL pw_axpy(v_rspace_embed(ispin), v_rspace_new(ispin), v_rspace_embed(ispin)%pw_grid%dvol)
1441 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1444 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1445 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1447 lri_v_int(ikind)%v_int = 0.0_dp
1449 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1450 lri_v_int, calculate_forces,
"LRI_AUX")
1452 CALL para_env%sum(lri_v_int(ikind)%v_int)
1454 IF (lri_env%exact_1c_terms)
THEN
1455 rho_ao => my_rho(ispin, :)
1456 ksmat => ks_matrix(ispin, :)
1457 CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), ksmat(1)%matrix, &
1458 rho_ao(1)%matrix, qs_env, &
1459 calculate_forces,
"ORB")
1461 IF (lri_env%ppl_ri)
THEN
1465 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1466 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1468 lri_v_int(ikind)%v_int = 0.0_dp
1470 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1471 lri_v_int, calculate_forces,
"RI_HXC")
1473 CALL para_env%sum(lri_v_int(ikind)%v_int)
1476 rho_ao => my_rho(ispin, :)
1477 ksmat => ks_matrix(ispin, :)
1478 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1479 pmat_kp=rho_ao, hmat_kp=ksmat, &
1481 calculate_forces=calculate_forces, &
1484 CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1487 SELECT CASE (dft_control%sic_method_id)
1490 CALL auxbas_pw_pool%give_back_pw(v_sic_rspace)
1491 DEALLOCATE (v_sic_rspace)
1493 DEALLOCATE (v_rspace_new)
1497 cpassert(dft_control%sic_method_id ==
sic_none)
1498 cpassert(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
1499 DO ispin = 1, nspins
1502 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1503 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace, dvol)
1506 IF (dft_control%do_sccs)
THEN
1507 CALL pw_axpy(v_sccs_rspace, v_rspace)
1510 IF (dft_control%apply_embed_pot)
THEN
1511 CALL pw_axpy(v_rspace_embed(ispin), v_rspace, v_rspace_embed(ispin)%pw_grid%dvol)
1512 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1515 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1516 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1518 lri_v_int(ikind)%v_int = 0.0_dp
1520 CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1521 lri_v_int, calculate_forces,
"LRI_AUX")
1523 CALL para_env%sum(lri_v_int(ikind)%v_int)
1525 IF (lri_env%exact_1c_terms)
THEN
1526 rho_ao => my_rho(ispin, :)
1527 ksmat => ks_matrix(ispin, :)
1528 CALL integrate_v_rspace_diagonal(v_rspace, ksmat(1)%matrix, &
1529 rho_ao(1)%matrix, qs_env, &
1530 calculate_forces,
"ORB")
1532 IF (lri_env%ppl_ri)
THEN
1536 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1537 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1539 lri_v_int(ikind)%v_int = 0.0_dp
1541 CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1542 lri_v_int, calculate_forces,
"RI_HXC")
1544 CALL para_env%sum(lri_v_int(ikind)%v_int)
1547 rho_ao => my_rho(ispin, :)
1548 ksmat => ks_matrix(ispin, :)
1549 CALL integrate_v_rspace(v_rspace=v_rspace, &
1553 calculate_forces=calculate_forces, &
1562 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1564 DO ispin = 1, nspins
1565 ksmat => ks_matrix(ispin, :)
1567 cell_to_index=cell_to_index)
1569 IF (calculate_forces)
THEN
1574 DO ispin = 1, nspins
1576 smat(1)%matrix, atomic_kind_set, ispin)
1578 IF (calculate_forces)
THEN
1579 rho_ao_nokp => rho_ao_kp(:, 1)
1584 IF (
ASSOCIATED(v_tau_rspace))
THEN
1585 IF (lrigpw .OR. rigpw)
THEN
1586 cpabort(
"LRIGPW/RIGPW not implemented for meta-GGAs")
1588 DO ispin = 1, nspins
1589 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1591 rho_ao => rho_ao_kp(ispin, :)
1592 ksmat => ks_matrix(ispin, :)
1593 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1594 pmat_kp=rho_ao, hmat_kp=ksmat, &
1596 calculate_forces=calculate_forces, compute_tau=.true., &
1598 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1600 DEALLOCATE (v_tau_rspace)
1604 IF (dft_control%do_admm)
THEN
1606 CALL get_admm_env(admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
1607 matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
1608 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
1609 IF (
ASSOCIATED(v_rspace_new_aux_fit))
THEN
1610 DO ispin = 1, nspins
1612 CALL pw_scale(v_rspace_new_aux_fit(ispin), v_rspace_new_aux_fit(ispin)%pw_grid%dvol)
1615 DO img = 1, dft_control%nimages
1616 CALL dbcsr_copy(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit(ispin, img)%matrix, &
1617 name=
"DFT exch. part of matrix_ks_aux_fit")
1625 CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
1626 basis_type =
"AUX_FIT"
1627 IF (admm_env%do_gapw)
THEN
1628 task_list => admm_env%admm_gapw_env%task_list
1629 basis_type =
"AUX_FIT_SOFT"
1633 IF (admm_env%do_admmp)
THEN
1634 fadm = admm_env%gsi(ispin)**2
1635 ELSE IF (admm_env%do_admms)
THEN
1636 fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
1639 rho_ao => rho_ao_aux(ispin, :)
1640 ksmat => matrix_ks_aux_fit(ispin, :)
1642 CALL integrate_v_rspace(v_rspace=v_rspace_new_aux_fit(ispin), &
1646 calculate_forces=calculate_forces, &
1649 basis_type=basis_type, &
1650 task_list_external=task_list)
1654 DO img = 1, dft_control%nimages
1655 CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
1656 matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
1659 CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
1661 DEALLOCATE (v_rspace_new_aux_fit)
1664 IF (
ASSOCIATED(v_tau_rspace_aux_fit))
THEN
1665 DO ispin = 1, nspins
1666 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace_aux_fit(ispin))
1668 DEALLOCATE (v_tau_rspace_aux_fit)
1672 IF (dft_control%apply_embed_pot)
DEALLOCATE (v_rspace_embed)
1674 CALL timestop(handle)
1696 REAL(kind=
dp) :: exc
1698 CHARACTER(*),
PARAMETER :: routinen =
'calculate_zmp_potential'
1700 INTEGER :: handle, my_val, nelectron, nspins
1701 INTEGER,
DIMENSION(2) :: nelectron_spin
1702 LOGICAL :: do_zmp_read, fermi_amaldi
1703 REAL(kind=
dp) :: lambda
1704 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tot_rho_ext_r
1717 CALL timeset(routinen, handle)
1718 NULLIFY (auxbas_pw_pool)
1720 NULLIFY (poisson_env)
1721 NULLIFY (v_rspace_new)
1722 NULLIFY (dft_control)
1723 NULLIFY (rho_r, rho_g, tot_rho_ext_r, rho_ext_g)
1729 nelectron_spin=nelectron_spin, &
1730 dft_control=dft_control)
1732 auxbas_pw_pool=auxbas_pw_pool, &
1733 poisson_env=poisson_env)
1734 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
1736 ALLOCATE (v_rspace_new(nspins))
1737 CALL auxbas_pw_pool%create_pw(pw=v_rspace_new(1))
1738 CALL auxbas_pw_pool%create_pw(pw=v_xc_rspace)
1741 do_zmp_read = dft_control%apply_external_vxc
1742 IF (do_zmp_read)
THEN
1743 CALL pw_copy(qs_env%external_vxc, v_rspace_new(1))
1745 v_rspace_new(1)%pw_grid%dvol
1748 REAL(kind=
dp) :: factor
1750 CALL auxbas_pw_pool%create_pw(pw=rho_eff_gspace)
1751 CALL auxbas_pw_pool%create_pw(pw=v_xc_gspace)
1758 tot_rho_r=tot_rho_ext_r)
1759 factor = tot_rho_ext_r(1)/factor
1761 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1762 CALL pw_axpy(rho_ext_g(1), rho_eff_gspace, alpha=-1.0_dp)
1768 CALL pw_scale(rho_eff_gspace, a=lambda)
1769 nelectron = nelectron_spin(1)
1770 factor = -1.0_dp/nelectron
1771 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1775 CALL pw_copy(v_rspace_new(1), v_xc_rspace)
1783 CALL auxbas_pw_pool%give_back_pw(rho_eff_gspace)
1784 CALL auxbas_pw_pool%give_back_pw(v_xc_gspace)
1788 CALL auxbas_pw_pool%give_back_pw(v_xc_rspace)
1790 CALL timestop(handle)
1805 TYPE(qs_environment_type),
POINTER :: qs_env
1806 TYPE(qs_rho_type),
POINTER :: rho
1807 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: v_rspace_embed
1808 TYPE(dft_control_type),
POINTER :: dft_control
1809 REAL(kind=dp) :: embed_corr
1810 LOGICAL :: just_energy
1812 CHARACTER(*),
PARAMETER :: routinen =
'get_embed_potential_energy'
1814 INTEGER :: handle, ispin
1815 REAL(kind=dp) :: embed_corr_local
1816 TYPE(pw_env_type),
POINTER :: pw_env
1817 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
1818 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r
1820 CALL timeset(routinen, handle)
1822 NULLIFY (auxbas_pw_pool)
1825 CALL get_qs_env(qs_env=qs_env, &
1828 CALL pw_env_get(pw_env=pw_env, &
1829 auxbas_pw_pool=auxbas_pw_pool)
1830 CALL qs_rho_get(rho, rho_r=rho_r)
1831 ALLOCATE (v_rspace_embed(dft_control%nspins))
1835 DO ispin = 1, dft_control%nspins
1836 CALL auxbas_pw_pool%create_pw(pw=v_rspace_embed(ispin))
1837 CALL pw_zero(v_rspace_embed(ispin))
1839 CALL pw_copy(qs_env%embed_pot, v_rspace_embed(ispin))
1840 embed_corr_local = 0.0_dp
1843 IF (dft_control%nspins .EQ. 2)
THEN
1844 IF (ispin .EQ. 1)
CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), 1.0_dp)
1845 IF (ispin .EQ. 2)
CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), -1.0_dp)
1848 embed_corr_local = pw_integral_ab(v_rspace_embed(ispin), rho_r(ispin))
1850 embed_corr = embed_corr + embed_corr_local
1855 IF (just_energy)
THEN
1856 DO ispin = 1, dft_control%nspins
1857 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1859 DEALLOCATE (v_rspace_embed)
1862 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.
Handles all functions related to the CELL.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
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_set(matrix, alpha)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
subroutine, public dbcsr_scale_by_vector(matrix, alpha, side)
Scales the rows/columns of given matrix.
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,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
subroutine, public cp_ddapc_apply_cd(qs_env, rho_tot_gspace, energy, v_hartree_gspace, calculate_forces, itype_of_density)
Routine to couple/decouple periodic images with the Bloechl scheme.
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_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 ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Utilities for hfx and admm methods.
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data, external_para_env)
Add the hfx contributions to the Hamiltonian.
Routines to calculate derivatives with respect to basis function origin.
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Types and set/get functions for HFX.
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Calculates integral matrices for LRIGPW method lri : local resolution of the identity.
subroutine, public v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
...
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Calculates forces for LRIGPW method lri : local resolution of the identity.
subroutine, public calculate_lri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the lri forces
subroutine, public calculate_ri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the ri forces
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
subroutine, public calculate_ri_ks_matrix(lri_env, lri_v_int, h_matrix, s_matrix, atomic_kind_set, ispin)
update of RIGPW KS matrix
Interface to the message passing library MPI.
Types containing essential information for running implicit (iterative) Poisson solver.
integer, parameter, public neumann_bc
integer, parameter, public mixed_bc
integer, parameter, public mixed_periodic_bc
integer, parameter, public periodic_bc
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
integer, parameter, public pw_poisson_implicit
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Defines CDFT control structures.
container for information about total charges on the grids
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
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.
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_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, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
subroutine, public qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale)
Modify the hartree potential in order to include the QM/MM correction.
subroutine, public get_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, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, 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, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public print_densities(qs_env, rho)
...
subroutine, public get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, just_energy)
...
subroutine, public low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, calculate_forces, auxbas_pw_pool)
do ROKS calculations yielding low spin states
subroutine, public sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, cdft_control, calculate_forces)
Sum up all potentials defined on the grid and integrate.
subroutine, public print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
Print detailed energies.
subroutine, public calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
Calculate the ZMP potential and energy as in Zhao, Morrison Parr PRA 50i, 2138 (1994) V_c^\lambda def...
subroutine, public calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on the spin density
subroutine, public sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on explicit orbitals
subroutine, public compute_matrix_vxc(qs_env, v_rspace, matrix_vxc)
compute matrix_vxc, defined via the potential created by qs_vxc_create ignores things like tau functi...
Definition and initialisation of the mo data type.
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.
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...
Exchange and Correlation functional calculations.
real(kind=dp) function, public xc_exc_calc(rho_r, rho_g, tau, xc_section, pw_pool)
calculates just the exchange and correlation energy (no vxc)
subroutine, public xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_section, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations. depending on the selected functional_routine calls ...
stores some data used in wavefunction fitting
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores some data used in construction of Kohn-Sham matrix
Contains information about kpoints.
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Container for information about total charges on the grids.
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.