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(:, :),
POINTER :: 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(:, :),
POINTER :: 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)
504 INTEGER,
INTENT(IN) :: nspins
505 INTEGER,
DIMENSION(2, 3),
INTENT(IN) :: bo
512 IF (needs%rho_1_3)
THEN
513 NULLIFY (rho_set%rho_1_3)
514 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
515 rho_set%owns%rho_1_3 = .true.
516 rho_set%has%rho_1_3 = .false.
520 NULLIFY (rho_set%rho)
521 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
522 rho_set%owns%rho = .true.
523 rho_set%has%rho = .false.
526 IF (needs%norm_drho)
THEN
527 NULLIFY (rho_set%norm_drho)
528 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
529 rho_set%owns%norm_drho = .true.
530 rho_set%has%norm_drho = .false.
535 NULLIFY (rho_set%drho(idir)%array)
536 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
538 rho_set%owns%drho = .true.
539 rho_set%has%drho = .false.
545 NULLIFY (rho_set%rho)
546 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
547 rho_set%owns%rho = .true.
548 rho_set%has%rho = .false.
551 IF (needs%rho_1_3)
THEN
552 NULLIFY (rho_set%rho_1_3)
553 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
554 rho_set%owns%rho_1_3 = .true.
555 rho_set%has%rho_1_3 = .false.
558 IF (needs%rho_spin_1_3)
THEN
559 NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
560 ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
561 ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
562 rho_set%owns%rho_spin_1_3 = .true.
563 rho_set%has%rho_spin_1_3 = .false.
566 IF (needs%rho_spin)
THEN
567 NULLIFY (rho_set%rhoa, rho_set%rhob)
568 ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
569 ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
570 rho_set%owns%rho_spin = .true.
571 rho_set%has%rho_spin = .false.
574 IF (needs%norm_drho)
THEN
575 NULLIFY (rho_set%norm_drho)
576 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
577 rho_set%owns%norm_drho = .true.
578 rho_set%has%norm_drho = .false.
581 IF (needs%norm_drho_spin)
THEN
582 NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
583 ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
584 ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
585 rho_set%owns%norm_drho_spin = .true.
586 rho_set%has%norm_drho_spin = .false.
591 NULLIFY (rho_set%drho(idir)%array)
592 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
594 rho_set%owns%drho = .true.
595 rho_set%has%drho = .false.
598 IF (needs%drho_spin)
THEN
600 NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
601 ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
602 ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
604 rho_set%owns%drho_spin = .true.
605 rho_set%has%drho_spin = .false.
612 NULLIFY (rho_set%tau)
613 ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
614 rho_set%owns%tau = .true.
616 IF (needs%tau_spin)
THEN
617 NULLIFY (rho_set%tau_a, rho_set%tau_b)
618 ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
619 ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
620 rho_set%owns%tau_spin = .true.
621 rho_set%has%tau_spin = .false.
625 IF (needs%laplace_rho)
THEN
626 NULLIFY (rho_set%laplace_rho)
627 ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
628 rho_set%owns%laplace_rho = .true.
630 IF (needs%laplace_rho_spin)
THEN
631 NULLIFY (rho_set%laplace_rhoa)
632 NULLIFY (rho_set%laplace_rhob)
633 ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
634 ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
635 rho_set%owns%laplace_rho_spin = .true.
636 rho_set%has%laplace_rho_spin = .true.
653 SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
656 LOGICAL,
INTENT(IN) :: lsd
657 INTEGER,
INTENT(IN) :: nspins
659 REAL(
dp),
DIMENSION(:, :, :),
INTENT(IN) :: rho
660 REAL(
dp),
DIMENSION(:, :, :, :),
INTENT(IN) :: drho
661 REAL(
dp),
DIMENSION(:, :, :),
INTENT(IN) :: tau
662 INTEGER,
INTENT(IN) :: na, ir
664 REAL(kind=
dp),
PARAMETER :: f13 = (1.0_dp/3.0_dp)
666 INTEGER :: ia, idir, my_nspins
667 LOGICAL :: gradient_f, tddft_split
670 tddft_split = .false.
671 IF (lsd .AND. nspins == 1)
THEN
679 cpassert(
SIZE(rho, 3) == 1)
681 SELECT CASE (my_nspins)
683 cpassert(.NOT. needs%rho_spin)
684 cpassert(.NOT. needs%drho_spin)
685 cpassert(.NOT. needs%norm_drho_spin)
686 cpassert(.NOT. needs%rho_spin_1_3)
689 cpabort(
"Unsupported number of spins")
692 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
693 needs%drho .OR. needs%norm_drho)
695 SELECT CASE (my_nspins)
698 IF (needs%rho_1_3)
THEN
700 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
702 rho_set%owns%rho_1_3 = .true.
703 rho_set%has%rho_1_3 = .true.
708 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
710 rho_set%owns%rho = .true.
711 rho_set%has%rho = .true.
714 IF (needs%norm_drho)
THEN
716 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
718 rho_set%owns%norm_drho = .true.
719 rho_set%has%norm_drho = .true.
725 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
728 rho_set%owns%drho = .true.
729 rho_set%has%drho = .true.
735 IF (.NOT. tddft_split)
THEN
737 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
741 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
744 rho_set%owns%rho = .true.
745 rho_set%has%rho = .true.
748 IF (needs%rho_1_3)
THEN
749 IF (.NOT. tddft_split)
THEN
751 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
755 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
758 rho_set%owns%rho_1_3 = .true.
759 rho_set%has%rho_1_3 = .true.
762 IF (needs%rho_spin_1_3)
THEN
763 IF (.NOT. tddft_split)
THEN
765 rho_set%rhoa_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
766 rho_set%rhob_1_3(ia, ir, 1) = max(rho(ia, ir, 2), 0.0_dp)**f13
770 rho_set%rhoa_1_3(ia, ir, 1) = max(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
771 rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
774 rho_set%owns%rho_spin_1_3 = .true.
775 rho_set%has%rho_spin_1_3 = .true.
778 IF (needs%rho_spin)
THEN
779 IF (.NOT. tddft_split)
THEN
781 rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
782 rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
786 rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
787 rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
790 rho_set%owns%rho_spin = .true.
791 rho_set%has%rho_spin = .true.
794 IF (needs%norm_drho)
THEN
795 IF (.NOT. tddft_split)
THEN
797 rho_set%norm_drho(ia, ir, 1) = sqrt( &
798 (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
799 (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
800 (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
804 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
807 rho_set%owns%norm_drho = .true.
808 rho_set%has%norm_drho = .true.
811 IF (needs%norm_drho_spin)
THEN
812 IF (.NOT. tddft_split)
THEN
814 rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
815 rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
819 rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
820 rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
823 rho_set%owns%norm_drho_spin = .true.
824 rho_set%has%norm_drho_spin = .true.
828 IF (.NOT. tddft_split)
THEN
831 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
837 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
841 rho_set%owns%drho = .true.
842 rho_set%has%drho = .true.
845 IF (needs%drho_spin)
THEN
846 IF (.NOT. tddft_split)
THEN
849 rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
850 rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
856 rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
857 rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
861 rho_set%owns%drho_spin = .true.
862 rho_set%has%drho_spin = .true.
868 IF (needs%tau .OR. needs%tau_spin)
THEN
869 cpassert(
SIZE(tau, 3) == my_nspins)
872 IF (my_nspins == 2)
THEN
874 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
876 rho_set%owns%tau = .true.
877 rho_set%has%tau = .true.
880 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
882 rho_set%owns%tau = .true.
883 rho_set%has%tau = .true.
886 IF (needs%tau_spin)
THEN
888 rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
889 rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
891 rho_set%owns%tau_spin = .true.
892 rho_set%has%tau_spin = .true.