61 SUBROUTINE vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, &
62 lsd, na, nr, exc, vxc, vxg, vtau, &
63 energy_only, adiabatic_rescale_factor)
78 INTEGER,
INTENT(in) :: deriv_order
80 REAL(
dp),
DIMENSION(:, :),
INTENT(IN) :: w
81 LOGICAL,
INTENT(IN) :: lsd
82 INTEGER,
INTENT(in) :: na, nr
84 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: vxc
85 REAL(
dp),
DIMENSION(:, :, :, :),
POINTER :: vxg
86 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: vtau
87 LOGICAL,
INTENT(IN),
OPTIONAL :: energy_only
88 REAL(
dp),
INTENT(IN),
OPTIONAL :: adiabatic_rescale_factor
90 CHARACTER(LEN=*),
PARAMETER :: routinen =
'vxc_of_r_new'
92 INTEGER :: handle, ia, idir, ir
93 LOGICAL :: gradient_f, my_only_energy
94 REAL(
dp) :: my_adiabatic_rescale_factor
95 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data
96 REAL(kind=
dp) :: drho_cutoff
99 CALL timeset(routinen, handle)
100 my_only_energy = .false.
101 IF (
PRESENT(energy_only)) my_only_energy = energy_only
103 IF (
PRESENT(adiabatic_rescale_factor))
THEN
104 my_adiabatic_rescale_factor = adiabatic_rescale_factor
106 my_adiabatic_rescale_factor = 1.0_dp
109 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
110 needs%drho .OR. needs%norm_drho)
116 deriv_set=deriv_set, &
117 deriv_order=deriv_order)
126 IF (
ASSOCIATED(deriv_att))
THEN
130 exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
136 IF (.NOT. my_only_energy)
THEN
140 IF (
ASSOCIATED(deriv_att))
THEN
142 vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
146 IF (
ASSOCIATED(deriv_att))
THEN
148 vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
152 IF (
ASSOCIATED(deriv_att))
THEN
154 vxc(:, :, 1) = vxc(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
155 vxc(:, :, 2) = vxc(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
160 IF (
ASSOCIATED(deriv_att))
THEN
162 vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
170 IF (
ASSOCIATED(deriv_att))
THEN
175 IF (rho_set%norm_drhoa(ia, ir, 1) > drho_cutoff)
THEN
176 vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
177 deriv_data(ia, ir, 1)*w(ia, ir)/ &
178 rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
180 vxg(idir, ia, ir, 1) = 0.0_dp
188 IF (
ASSOCIATED(deriv_att))
THEN
193 IF (rho_set%norm_drhob(ia, ir, 1) > drho_cutoff)
THEN
194 vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
195 deriv_data(ia, ir, 1)*w(ia, ir)/ &
196 rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
198 vxg(idir, ia, ir, 2) = 0.0_dp
207 IF (
ASSOCIATED(deriv_att))
THEN
212 IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff)
THEN
213 vxg(idir, ia, ir, 1:2) = &
214 vxg(idir, ia, ir, 1:2) + ( &
215 rho_set%drhoa(idir)%array(ia, ir, 1) + &
216 rho_set%drhob(idir)%array(ia, ir, 1))* &
217 deriv_data(ia, ir, 1)*w(ia, ir)/rho_set%norm_drho(ia, ir, 1)* &
218 my_adiabatic_rescale_factor
227 IF (
ASSOCIATED(deriv_att))
THEN
231 IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff)
THEN
233 vxg(idir, ia, ir, 1) = rho_set%drho(idir)%array(ia, ir, 1)* &
234 deriv_data(ia, ir, 1)*w(ia, ir)/ &
235 rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
238 vxg(1:3, ia, ir, 1) = 0.0_dp
248 IF (
ASSOCIATED(deriv_att))
THEN
250 vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
254 IF (
ASSOCIATED(deriv_att))
THEN
256 vtau(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
260 IF (
ASSOCIATED(deriv_att))
THEN
262 vtau(:, :, 1) = vtau(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
263 vtau(:, :, 2) = vtau(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
268 IF (
ASSOCIATED(deriv_att))
THEN
270 vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
276 CALL timestop(handle)
295 SUBROUTINE vxc_of_r_epr(xc_fun_section, rho_set, deriv_set, needs, w, &
296 lsd, na, nr, exc, vxc, vxg, vtau)
302 REAL(
dp),
DIMENSION(:, :),
INTENT(IN) :: w
303 LOGICAL,
INTENT(IN) :: lsd
304 INTEGER,
INTENT(in) :: na, nr
306 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: vxc
307 REAL(
dp),
DIMENSION(:, :, :, :),
POINTER :: vxg
308 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: vtau
310 CHARACTER(LEN=*),
PARAMETER :: routinen =
'vxc_of_r_epr'
312 INTEGER :: handle, ia, idir, ir, my_deriv_order
313 LOGICAL :: gradient_f
314 REAL(
dp) :: my_adiabatic_rescale_factor
315 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data
316 REAL(kind=
dp) :: drho_cutoff
319 CALL timeset(routinen, handle)
324 my_adiabatic_rescale_factor = 1.0_dp
327 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
328 needs%drho .OR. needs%norm_drho)
334 deriv_set=deriv_set, &
335 deriv_order=my_deriv_order)
345 IF (
ASSOCIATED(deriv_att))
THEN
350 vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
351 deriv_data(ia, ir, 1)
358 IF (
ASSOCIATED(deriv_att))
THEN
363 vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
364 deriv_data(ia, ir, 1)
374 IF (
ASSOCIATED(deriv_att))
THEN
378 exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
384 CALL timestop(handle)
517 INTEGER,
INTENT(IN) :: nspins
518 INTEGER,
DIMENSION(2, 3),
INTENT(IN) :: bo
525 IF (needs%rho_1_3)
THEN
526 NULLIFY (rho_set%rho_1_3)
527 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
528 rho_set%owns%rho_1_3 = .true.
529 rho_set%has%rho_1_3 = .false.
533 NULLIFY (rho_set%rho)
534 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
535 rho_set%owns%rho = .true.
536 rho_set%has%rho = .false.
539 IF (needs%norm_drho)
THEN
540 NULLIFY (rho_set%norm_drho)
541 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
542 rho_set%owns%norm_drho = .true.
543 rho_set%has%norm_drho = .false.
548 NULLIFY (rho_set%drho(idir)%array)
549 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
551 rho_set%owns%drho = .true.
552 rho_set%has%drho = .false.
558 NULLIFY (rho_set%rho)
559 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
560 rho_set%owns%rho = .true.
561 rho_set%has%rho = .false.
564 IF (needs%rho_1_3)
THEN
565 NULLIFY (rho_set%rho_1_3)
566 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
567 rho_set%owns%rho_1_3 = .true.
568 rho_set%has%rho_1_3 = .false.
571 IF (needs%rho_spin_1_3)
THEN
572 NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
573 ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
574 ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
575 rho_set%owns%rho_spin_1_3 = .true.
576 rho_set%has%rho_spin_1_3 = .false.
579 IF (needs%rho_spin)
THEN
580 NULLIFY (rho_set%rhoa, rho_set%rhob)
581 ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
582 ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
583 rho_set%owns%rho_spin = .true.
584 rho_set%has%rho_spin = .false.
587 IF (needs%norm_drho)
THEN
588 NULLIFY (rho_set%norm_drho)
589 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
590 rho_set%owns%norm_drho = .true.
591 rho_set%has%norm_drho = .false.
594 IF (needs%norm_drho_spin)
THEN
595 NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
596 ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
597 ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
598 rho_set%owns%norm_drho_spin = .true.
599 rho_set%has%norm_drho_spin = .false.
604 NULLIFY (rho_set%drho(idir)%array)
605 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
607 rho_set%owns%drho = .true.
608 rho_set%has%drho = .false.
611 IF (needs%drho_spin)
THEN
613 NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
614 ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
615 ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
617 rho_set%owns%drho_spin = .true.
618 rho_set%has%drho_spin = .false.
625 NULLIFY (rho_set%tau)
626 ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
627 rho_set%owns%tau = .true.
629 IF (needs%tau_spin)
THEN
630 NULLIFY (rho_set%tau_a, rho_set%tau_b)
631 ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
632 ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
633 rho_set%owns%tau_spin = .true.
634 rho_set%has%tau_spin = .false.
638 IF (needs%laplace_rho)
THEN
639 NULLIFY (rho_set%laplace_rho)
640 ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
641 rho_set%owns%laplace_rho = .true.
643 IF (needs%laplace_rho_spin)
THEN
644 NULLIFY (rho_set%laplace_rhoa)
645 NULLIFY (rho_set%laplace_rhob)
646 ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
647 ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
648 rho_set%owns%laplace_rho_spin = .true.
649 rho_set%has%laplace_rho_spin = .true.
666 SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
669 LOGICAL,
INTENT(IN) :: lsd
670 INTEGER,
INTENT(IN) :: nspins
672 REAL(
dp),
DIMENSION(:, :, :),
INTENT(IN) :: rho
673 REAL(
dp),
DIMENSION(:, :, :, :),
INTENT(IN) :: drho
674 REAL(
dp),
DIMENSION(:, :, :),
INTENT(IN) :: tau
675 INTEGER,
INTENT(IN) :: na, ir
677 REAL(kind=
dp),
PARAMETER :: f13 = (1.0_dp/3.0_dp)
679 INTEGER :: ia, idir, my_nspins
680 LOGICAL :: gradient_f, tddft_split
683 tddft_split = .false.
684 IF (lsd .AND. nspins == 1)
THEN
692 cpassert(
SIZE(rho, 3) == 1)
694 SELECT CASE (my_nspins)
696 cpassert(.NOT. needs%rho_spin)
697 cpassert(.NOT. needs%drho_spin)
698 cpassert(.NOT. needs%norm_drho_spin)
699 cpassert(.NOT. needs%rho_spin_1_3)
702 cpabort(
"Unsupported number of spins")
705 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
706 needs%drho .OR. needs%norm_drho)
708 SELECT CASE (my_nspins)
711 IF (needs%rho_1_3)
THEN
713 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
715 rho_set%owns%rho_1_3 = .true.
716 rho_set%has%rho_1_3 = .true.
721 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
723 rho_set%owns%rho = .true.
724 rho_set%has%rho = .true.
727 IF (needs%norm_drho)
THEN
729 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
731 rho_set%owns%norm_drho = .true.
732 rho_set%has%norm_drho = .true.
738 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
741 rho_set%owns%drho = .true.
742 rho_set%has%drho = .true.
748 IF (.NOT. tddft_split)
THEN
750 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
754 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
757 rho_set%owns%rho = .true.
758 rho_set%has%rho = .true.
761 IF (needs%rho_1_3)
THEN
762 IF (.NOT. tddft_split)
THEN
764 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
768 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
771 rho_set%owns%rho_1_3 = .true.
772 rho_set%has%rho_1_3 = .true.
775 IF (needs%rho_spin_1_3)
THEN
776 IF (.NOT. tddft_split)
THEN
778 rho_set%rhoa_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
779 rho_set%rhob_1_3(ia, ir, 1) = max(rho(ia, ir, 2), 0.0_dp)**f13
783 rho_set%rhoa_1_3(ia, ir, 1) = max(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
784 rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
787 rho_set%owns%rho_spin_1_3 = .true.
788 rho_set%has%rho_spin_1_3 = .true.
791 IF (needs%rho_spin)
THEN
792 IF (.NOT. tddft_split)
THEN
794 rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
795 rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
799 rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
800 rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
803 rho_set%owns%rho_spin = .true.
804 rho_set%has%rho_spin = .true.
807 IF (needs%norm_drho)
THEN
808 IF (.NOT. tddft_split)
THEN
810 rho_set%norm_drho(ia, ir, 1) = sqrt( &
811 (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
812 (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
813 (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
817 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
820 rho_set%owns%norm_drho = .true.
821 rho_set%has%norm_drho = .true.
824 IF (needs%norm_drho_spin)
THEN
825 IF (.NOT. tddft_split)
THEN
827 rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
828 rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
832 rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
833 rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
836 rho_set%owns%norm_drho_spin = .true.
837 rho_set%has%norm_drho_spin = .true.
841 IF (.NOT. tddft_split)
THEN
844 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
850 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
854 rho_set%owns%drho = .true.
855 rho_set%has%drho = .true.
858 IF (needs%drho_spin)
THEN
859 IF (.NOT. tddft_split)
THEN
862 rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
863 rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
869 rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
870 rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
874 rho_set%owns%drho_spin = .true.
875 rho_set%has%drho_spin = .true.
881 IF (needs%tau .OR. needs%tau_spin)
THEN
882 cpassert(
SIZE(tau, 3) == my_nspins)
885 IF (my_nspins == 2)
THEN
887 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
889 rho_set%owns%tau = .true.
890 rho_set%has%tau = .true.
893 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
895 rho_set%owns%tau = .true.
896 rho_set%has%tau = .true.
899 IF (needs%tau_spin)
THEN
901 rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
902 rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
904 rho_set%owns%tau_spin = .true.
905 rho_set%has%tau_spin = .true.