62 SUBROUTINE vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, &
63 lsd, na, nr, exc, vxc, vxg, vtau, &
64 energy_only, epr_xc, adiabatic_rescale_factor)
79 INTEGER,
INTENT(in) :: deriv_order
81 REAL(
dp),
DIMENSION(:, :),
POINTER :: w
82 LOGICAL,
INTENT(IN) :: lsd
83 INTEGER,
INTENT(in) :: na, nr
85 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: vxc
86 REAL(
dp),
DIMENSION(:, :, :, :),
POINTER :: vxg
87 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: vtau
88 LOGICAL,
INTENT(IN),
OPTIONAL :: energy_only, epr_xc
89 REAL(
dp),
INTENT(IN),
OPTIONAL :: adiabatic_rescale_factor
91 CHARACTER(LEN=*),
PARAMETER :: routinen =
'vxc_of_r_new'
93 INTEGER :: handle, ia, idir, ir, my_deriv_order
94 LOGICAL :: gradient_f, my_epr_xc, my_only_energy
95 REAL(
dp) :: my_adiabatic_rescale_factor
96 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: deriv_data
97 REAL(kind=
dp) :: drho_cutoff
100 CALL timeset(routinen, handle)
101 my_only_energy = .false.
102 IF (
PRESENT(energy_only)) my_only_energy = energy_only
104 IF (
PRESENT(adiabatic_rescale_factor))
THEN
105 my_adiabatic_rescale_factor = adiabatic_rescale_factor
107 my_adiabatic_rescale_factor = 1.0_dp
112 IF (
PRESENT(epr_xc)) my_epr_xc = epr_xc
113 my_deriv_order = deriv_order
114 IF (my_epr_xc) my_deriv_order = 2
116 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
117 needs%drho .OR. needs%norm_drho)
123 deriv_set=deriv_set, &
124 deriv_order=my_deriv_order)
135 IF (
ASSOCIATED(deriv_att))
THEN
140 vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
141 deriv_data(ia, ir, 1)
148 IF (
ASSOCIATED(deriv_att))
THEN
153 vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
154 deriv_data(ia, ir, 1)
164 IF (
ASSOCIATED(deriv_att))
THEN
168 exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
177 IF (
ASSOCIATED(deriv_att))
THEN
181 exc = exc + deriv_data(ia, ir, 1)*w(ia, ir)
187 IF (.NOT. my_only_energy)
THEN
191 IF (
ASSOCIATED(deriv_att))
THEN
193 vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
197 IF (
ASSOCIATED(deriv_att))
THEN
199 vxc(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
203 IF (
ASSOCIATED(deriv_att))
THEN
205 vxc(:, :, 1) = vxc(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
206 vxc(:, :, 2) = vxc(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
211 IF (
ASSOCIATED(deriv_att))
THEN
213 vxc(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
221 IF (
ASSOCIATED(deriv_att))
THEN
226 IF (rho_set%norm_drhoa(ia, ir, 1) > drho_cutoff)
THEN
227 vxg(idir, ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)* &
228 deriv_data(ia, ir, 1)*w(ia, ir)/ &
229 rho_set%norm_drhoa(ia, ir, 1)*my_adiabatic_rescale_factor
231 vxg(idir, ia, ir, 1) = 0.0_dp
239 IF (
ASSOCIATED(deriv_att))
THEN
244 IF (rho_set%norm_drhob(ia, ir, 1) > drho_cutoff)
THEN
245 vxg(idir, ia, ir, 2) = rho_set%drhob(idir)%array(ia, ir, 1)* &
246 deriv_data(ia, ir, 1)*w(ia, ir)/ &
247 rho_set%norm_drhob(ia, ir, 1)*my_adiabatic_rescale_factor
249 vxg(idir, ia, ir, 2) = 0.0_dp
258 IF (
ASSOCIATED(deriv_att))
THEN
263 IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff)
THEN
264 vxg(idir, ia, ir, 1:2) = &
265 vxg(idir, ia, ir, 1:2) + ( &
266 rho_set%drhoa(idir)%array(ia, ir, 1) + &
267 rho_set%drhob(idir)%array(ia, ir, 1))* &
268 deriv_data(ia, ir, 1)*w(ia, ir)/rho_set%norm_drho(ia, ir, 1)* &
269 my_adiabatic_rescale_factor
278 IF (
ASSOCIATED(deriv_att))
THEN
282 IF (rho_set%norm_drho(ia, ir, 1) > drho_cutoff)
THEN
284 vxg(idir, ia, ir, 1) = rho_set%drho(idir)%array(ia, ir, 1)* &
285 deriv_data(ia, ir, 1)*w(ia, ir)/ &
286 rho_set%norm_drho(ia, ir, 1)*my_adiabatic_rescale_factor
289 vxg(1:3, ia, ir, 1) = 0.0_dp
299 IF (
ASSOCIATED(deriv_att))
THEN
301 vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
305 IF (
ASSOCIATED(deriv_att))
THEN
307 vtau(:, :, 2) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
311 IF (
ASSOCIATED(deriv_att))
THEN
313 vtau(:, :, 1) = vtau(:, :, 1) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
314 vtau(:, :, 2) = vtau(:, :, 2) + deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
319 IF (
ASSOCIATED(deriv_att))
THEN
321 vtau(:, :, 1) = deriv_data(:, :, 1)*w(:, :)*my_adiabatic_rescale_factor
328 CALL timestop(handle)
442 INTEGER,
INTENT(IN) :: nspins
443 INTEGER,
DIMENSION(2, 3),
INTENT(IN) :: bo
450 IF (needs%rho_1_3)
THEN
451 NULLIFY (rho_set%rho_1_3)
452 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
453 rho_set%owns%rho_1_3 = .true.
454 rho_set%has%rho_1_3 = .false.
458 NULLIFY (rho_set%rho)
459 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
460 rho_set%owns%rho = .true.
461 rho_set%has%rho = .false.
464 IF (needs%norm_drho)
THEN
465 NULLIFY (rho_set%norm_drho)
466 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
467 rho_set%owns%norm_drho = .true.
468 rho_set%has%norm_drho = .false.
473 NULLIFY (rho_set%drho(idir)%array)
474 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
476 rho_set%owns%drho = .true.
477 rho_set%has%drho = .false.
483 NULLIFY (rho_set%rho)
484 ALLOCATE (rho_set%rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
485 rho_set%owns%rho = .true.
486 rho_set%has%rho = .false.
489 IF (needs%rho_1_3)
THEN
490 NULLIFY (rho_set%rho_1_3)
491 ALLOCATE (rho_set%rho_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
492 rho_set%owns%rho_1_3 = .true.
493 rho_set%has%rho_1_3 = .false.
496 IF (needs%rho_spin_1_3)
THEN
497 NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
498 ALLOCATE (rho_set%rhoa_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
499 ALLOCATE (rho_set%rhob_1_3(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
500 rho_set%owns%rho_spin_1_3 = .true.
501 rho_set%has%rho_spin_1_3 = .false.
504 IF (needs%rho_spin)
THEN
505 NULLIFY (rho_set%rhoa, rho_set%rhob)
506 ALLOCATE (rho_set%rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
507 ALLOCATE (rho_set%rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
508 rho_set%owns%rho_spin = .true.
509 rho_set%has%rho_spin = .false.
512 IF (needs%norm_drho)
THEN
513 NULLIFY (rho_set%norm_drho)
514 ALLOCATE (rho_set%norm_drho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
515 rho_set%owns%norm_drho = .true.
516 rho_set%has%norm_drho = .false.
519 IF (needs%norm_drho_spin)
THEN
520 NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
521 ALLOCATE (rho_set%norm_drhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
522 ALLOCATE (rho_set%norm_drhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
523 rho_set%owns%norm_drho_spin = .true.
524 rho_set%has%norm_drho_spin = .false.
529 NULLIFY (rho_set%drho(idir)%array)
530 ALLOCATE (rho_set%drho(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
532 rho_set%owns%drho = .true.
533 rho_set%has%drho = .false.
536 IF (needs%drho_spin)
THEN
538 NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
539 ALLOCATE (rho_set%drhoa(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
540 ALLOCATE (rho_set%drhob(idir)%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
542 rho_set%owns%drho_spin = .true.
543 rho_set%has%drho_spin = .false.
550 NULLIFY (rho_set%tau)
551 ALLOCATE (rho_set%tau(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
552 rho_set%owns%tau = .true.
554 IF (needs%tau_spin)
THEN
555 NULLIFY (rho_set%tau_a, rho_set%tau_b)
556 ALLOCATE (rho_set%tau_a(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
557 ALLOCATE (rho_set%tau_b(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
558 rho_set%owns%tau_spin = .true.
559 rho_set%has%tau_spin = .false.
563 IF (needs%laplace_rho)
THEN
564 NULLIFY (rho_set%laplace_rho)
565 ALLOCATE (rho_set%laplace_rho(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
566 rho_set%owns%laplace_rho = .true.
568 IF (needs%laplace_rho_spin)
THEN
569 NULLIFY (rho_set%laplace_rhoa)
570 NULLIFY (rho_set%laplace_rhob)
571 ALLOCATE (rho_set%laplace_rhoa(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
572 ALLOCATE (rho_set%laplace_rhob(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
573 rho_set%owns%laplace_rho_spin = .true.
574 rho_set%has%laplace_rho_spin = .true.
591 SUBROUTINE fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
594 LOGICAL,
INTENT(IN) :: lsd
595 INTEGER,
INTENT(IN) :: nspins
597 REAL(
dp),
DIMENSION(:, :, :),
INTENT(IN) :: rho
598 REAL(
dp),
DIMENSION(:, :, :, :),
INTENT(IN) :: drho
599 REAL(
dp),
DIMENSION(:, :, :),
INTENT(IN) :: tau
600 INTEGER,
INTENT(IN) :: na, ir
602 REAL(kind=
dp),
PARAMETER :: f13 = (1.0_dp/3.0_dp)
604 INTEGER :: ia, idir, my_nspins
605 LOGICAL :: gradient_f, tddft_split
608 tddft_split = .false.
609 IF (lsd .AND. nspins == 1)
THEN
617 cpassert(
SIZE(rho, 3) == 1)
619 SELECT CASE (my_nspins)
621 cpassert(.NOT. needs%rho_spin)
622 cpassert(.NOT. needs%drho_spin)
623 cpassert(.NOT. needs%norm_drho_spin)
624 cpassert(.NOT. needs%rho_spin_1_3)
627 cpabort(
"Unsupported number of spins")
630 gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
631 needs%drho .OR. needs%norm_drho)
633 SELECT CASE (my_nspins)
636 IF (needs%rho_1_3)
THEN
638 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
640 rho_set%owns%rho_1_3 = .true.
641 rho_set%has%rho_1_3 = .true.
646 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
648 rho_set%owns%rho = .true.
649 rho_set%has%rho = .true.
652 IF (needs%norm_drho)
THEN
654 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
656 rho_set%owns%norm_drho = .true.
657 rho_set%has%norm_drho = .true.
663 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
666 rho_set%owns%drho = .true.
667 rho_set%has%drho = .true.
673 IF (.NOT. tddft_split)
THEN
675 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1) + rho(ia, ir, 2)
679 rho_set%rho(ia, ir, 1) = rho(ia, ir, 1)
682 rho_set%owns%rho = .true.
683 rho_set%has%rho = .true.
686 IF (needs%rho_1_3)
THEN
687 IF (.NOT. tddft_split)
THEN
689 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1) + rho(ia, ir, 2), 0.0_dp)**f13
693 rho_set%rho_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
696 rho_set%owns%rho_1_3 = .true.
697 rho_set%has%rho_1_3 = .true.
700 IF (needs%rho_spin_1_3)
THEN
701 IF (.NOT. tddft_split)
THEN
703 rho_set%rhoa_1_3(ia, ir, 1) = max(rho(ia, ir, 1), 0.0_dp)**f13
704 rho_set%rhob_1_3(ia, ir, 1) = max(rho(ia, ir, 2), 0.0_dp)**f13
708 rho_set%rhoa_1_3(ia, ir, 1) = max(0.5_dp*rho(ia, ir, 1), 0.0_dp)**f13
709 rho_set%rhob_1_3(ia, ir, 1) = rho_set%rhoa_1_3(ia, ir, 1)
712 rho_set%owns%rho_spin_1_3 = .true.
713 rho_set%has%rho_spin_1_3 = .true.
716 IF (needs%rho_spin)
THEN
717 IF (.NOT. tddft_split)
THEN
719 rho_set%rhoa(ia, ir, 1) = rho(ia, ir, 1)
720 rho_set%rhob(ia, ir, 1) = rho(ia, ir, 2)
724 rho_set%rhoa(ia, ir, 1) = 0.5_dp*rho(ia, ir, 1)
725 rho_set%rhob(ia, ir, 1) = rho_set%rhoa(ia, ir, 1)
728 rho_set%owns%rho_spin = .true.
729 rho_set%has%rho_spin = .true.
732 IF (needs%norm_drho)
THEN
733 IF (.NOT. tddft_split)
THEN
735 rho_set%norm_drho(ia, ir, 1) = sqrt( &
736 (drho(1, ia, ir, 1) + drho(1, ia, ir, 2))**2 + &
737 (drho(2, ia, ir, 1) + drho(2, ia, ir, 2))**2 + &
738 (drho(3, ia, ir, 1) + drho(3, ia, ir, 2))**2)
742 rho_set%norm_drho(ia, ir, 1) = drho(4, ia, ir, 1)
745 rho_set%owns%norm_drho = .true.
746 rho_set%has%norm_drho = .true.
749 IF (needs%norm_drho_spin)
THEN
750 IF (.NOT. tddft_split)
THEN
752 rho_set%norm_drhoa(ia, ir, 1) = drho(4, ia, ir, 1)
753 rho_set%norm_drhob(ia, ir, 1) = drho(4, ia, ir, 2)
757 rho_set%norm_drhoa(ia, ir, 1) = 0.5_dp*drho(4, ia, ir, 1)
758 rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1)
761 rho_set%owns%norm_drho_spin = .true.
762 rho_set%has%norm_drho_spin = .true.
766 IF (.NOT. tddft_split)
THEN
769 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1) + drho(idir, ia, ir, 2)
775 rho_set%drho(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
779 rho_set%owns%drho = .true.
780 rho_set%has%drho = .true.
783 IF (needs%drho_spin)
THEN
784 IF (.NOT. tddft_split)
THEN
787 rho_set%drhoa(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 1)
788 rho_set%drhob(idir)%array(ia, ir, 1) = drho(idir, ia, ir, 2)
794 rho_set%drhoa(idir)%array(ia, ir, 1) = 0.5_dp*drho(idir, ia, ir, 1)
795 rho_set%drhob(idir)%array(ia, ir, 1) = rho_set%drhoa(idir)%array(ia, ir, 1)
799 rho_set%owns%drho_spin = .true.
800 rho_set%has%drho_spin = .true.
806 IF (needs%tau .OR. needs%tau_spin)
THEN
807 cpassert(
SIZE(tau, 3) == my_nspins)
810 IF (my_nspins == 2)
THEN
812 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1) + tau(ia, ir, 2)
814 rho_set%owns%tau = .true.
815 rho_set%has%tau = .true.
818 rho_set%tau(ia, ir, 1) = tau(ia, ir, 1)
820 rho_set%owns%tau = .true.
821 rho_set%has%tau = .true.
824 IF (needs%tau_spin)
THEN
826 rho_set%tau_a(ia, ir, 1) = tau(ia, ir, 1)
827 rho_set%tau_b(ia, ir, 1) = tau(ia, ir, 2)
829 rho_set%owns%tau_spin = .true.
830 rho_set%has%tau_spin = .true.