113#include "./base/base_uses.f90"
119 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_collocate_density'
139 MODULE PROCEDURE calculate_rho_core_r3d_rs
140 MODULE PROCEDURE calculate_rho_core_c1d_gs
144 MODULE PROCEDURE calculate_rho_resp_all_r3d_rs, calculate_rho_resp_all_c1d_gs
159 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_nlcc'
161 INTEGER :: atom_a, handle, iatom, iexp_nlcc, ikind, &
162 ithread, j, n, natom, nc, nexp_nlcc, &
163 ni, npme, nthread, subpatch_pattern
164 INTEGER,
DIMENSION(:),
POINTER :: atom_list, cores, nct_nlcc
166 REAL(kind=
dp) :: alpha, eps_rho_rspace, radius
167 REAL(kind=
dp),
DIMENSION(3) :: ra
168 REAL(kind=
dp),
DIMENSION(:),
POINTER :: alpha_nlcc
169 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: cval_nlcc, pab
177 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
180 CALL timeset(routinen, handle)
182 NULLIFY (cell, dft_control, pab, particle_set, atomic_kind_set, &
183 qs_kind_set, atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
186 atomic_kind_set=atomic_kind_set, &
187 qs_kind_set=qs_kind_set, &
189 dft_control=dft_control, &
190 particle_set=particle_set, &
192 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
193 auxbas_pw_pool=auxbas_pw_pool)
197 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
199 DO ikind = 1,
SIZE(atomic_kind_set)
200 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
201 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
203 IF (.NOT.
ASSOCIATED(gth_potential)) cycle
204 CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
205 alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
207 IF (.NOT. nlcc) cycle
209 DO iexp_nlcc = 1, nexp_nlcc
211 alpha = alpha_nlcc(iexp_nlcc)
212 nc = nct_nlcc(iexp_nlcc)
215 ALLOCATE (pab(ni, 1))
229 pab(1, 1) = cval_nlcc(1, iexp_nlcc)
232 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
234 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
236 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
239 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
241 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
243 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
245 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
247 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
249 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
252 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
254 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
256 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
258 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
260 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
262 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
264 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
266 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
268 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
270 pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
275 IF (dft_control%nspins == 2) pab = pab*0.5_dp
278 atom_a = atom_list(iatom)
279 ra(:) =
pbc(particle_set(atom_a)%r, cell)
280 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
282 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
295 atom_a = atom_list(iatom)
296 ra(:) =
pbc(particle_set(atom_a)%r, cell)
300 ra=ra, rb=ra, rp=ra, &
301 zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
302 pab=pab, o1=0, o2=0, &
303 prefactor=1.0_dp, cutoff=0.0_dp)
306 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
308 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
318 IF (
ASSOCIATED(cores))
THEN
324 CALL timestop(handle)
338 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_ppl_grid'
340 INTEGER :: atom_a, handle, iatom, ikind, ithread, &
341 j, lppl, n, natom, ni, npme, nthread, &
343 INTEGER,
DIMENSION(:),
POINTER :: atom_list, cores
344 REAL(kind=
dp) :: alpha, eps_rho_rspace, radius
345 REAL(kind=
dp),
DIMENSION(3) :: ra
346 REAL(kind=
dp),
DIMENSION(:),
POINTER :: cexp_ppl
347 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
355 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
358 CALL timeset(routinen, handle)
360 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
361 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
364 atomic_kind_set=atomic_kind_set, &
365 qs_kind_set=qs_kind_set, &
367 dft_control=dft_control, &
368 particle_set=particle_set, &
370 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
371 auxbas_pw_pool=auxbas_pw_pool)
375 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
377 DO ikind = 1,
SIZE(atomic_kind_set)
378 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
379 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
381 IF (.NOT.
ASSOCIATED(gth_potential)) cycle
382 CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
387 ALLOCATE (pab(ni, 1))
401 pab(1, 1) = cexp_ppl(1)
404 pab(n, 1) = cexp_ppl(2)
406 pab(n, 1) = cexp_ppl(2)
408 pab(n, 1) = cexp_ppl(2)
411 pab(n, 1) = cexp_ppl(3)
413 pab(n, 1) = cexp_ppl(3)
415 pab(n, 1) = cexp_ppl(3)
417 pab(n, 1) = 2._dp*cexp_ppl(3)
419 pab(n, 1) = 2._dp*cexp_ppl(3)
421 pab(n, 1) = 2._dp*cexp_ppl(3)
424 pab(n, 1) = cexp_ppl(4)
426 pab(n, 1) = cexp_ppl(4)
428 pab(n, 1) = cexp_ppl(4)
430 pab(n, 1) = 3._dp*cexp_ppl(4)
432 pab(n, 1) = 3._dp*cexp_ppl(4)
434 pab(n, 1) = 3._dp*cexp_ppl(4)
436 pab(n, 1) = 3._dp*cexp_ppl(4)
438 pab(n, 1) = 3._dp*cexp_ppl(4)
440 pab(n, 1) = 3._dp*cexp_ppl(4)
442 pab(n, 1) = 6._dp*cexp_ppl(4)
449 atom_a = atom_list(iatom)
450 ra(:) =
pbc(particle_set(atom_a)%r, cell)
451 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
453 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
463 IF (npme .GT. 0)
THEN
467 atom_a = atom_list(iatom)
468 ra(:) =
pbc(particle_set(atom_a)%r, cell)
473 lb_min=0, lb_max=0, &
474 ra=ra, rb=ra, rp=ra, &
475 zetp=alpha, eps=eps_rho_rspace, &
476 pab=pab, o1=0, o2=0, &
477 prefactor=1.0_dp, cutoff=0.0_dp)
480 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
482 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
491 IF (
ASSOCIATED(cores))
THEN
497 CALL timestop(handle)
517 lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
523 REAL(kind=
dp),
INTENT(OUT) :: total_rho
524 CHARACTER(len=*),
INTENT(IN) :: basis_type
525 LOGICAL,
INTENT(IN) :: exact_1c_terms
527 INTEGER,
DIMENSION(:),
OPTIONAL :: atomlist
529 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_lri_rho_elec'
531 INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
532 m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, nseta, offset, sgfa, sgfb
533 INTEGER,
DIMENSION(:),
POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
534 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa
536 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: map_it
537 LOGICAL,
ALLOCATABLE,
DIMENSION(:, :) :: map_it2
538 REAL(kind=
dp) :: eps_rho_rspace, radius, zetp
539 REAL(kind=
dp),
DIMENSION(3) :: ra
540 REAL(kind=
dp),
DIMENSION(:),
POINTER :: aci
541 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: p_block, pab, sphi_a, work, zeta
552 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
556 NULLIFY (aci, atomic_kind_set, qs_kind_set, atom_list, cell, &
557 dft_control, first_sgfa, gridlevel_info, la_max, &
558 la_min, lri_basis_set, npgfa, nsgfa, &
559 pab, particle_set, pw_env, pw_pools, rs_grid, rs_rho, sphi_a, &
562 CALL timeset(routinen, handle)
564 IF (exact_1c_terms)
THEN
565 cpassert(
PRESENT(pmat))
568 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
569 atomic_kind_set=atomic_kind_set, &
570 cell=cell, particle_set=particle_set, &
572 dft_control=dft_control)
574 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
575 gridlevel_info => pw_env%gridlevel_info
578 cpassert(
ASSOCIATED(pw_env))
579 CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, pw_pools=pw_pools)
586 DO igrid_level = 1, gridlevel_info%ngrid_levels
592 maxco=maxco, basis_type=basis_type)
594 ALLOCATE (pab(maxco, 1))
596 my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
597 group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
599 DO ikind = 1,
SIZE(atomic_kind_set)
601 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
602 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
606 lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
607 sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
610 atom_a = atom_list(iatom)
611 IF (
PRESENT(atomlist))
THEN
612 IF (atomlist(atom_a) == 0) cycle
614 ra(:) =
pbc(particle_set(atom_a)%r, cell)
615 aci => lri_coef(ikind)%acoef(iatom, :)
617 m1 = maxval(npgfa(1:nseta))
618 ALLOCATE (map_it(m1))
622 DO ipgf = 1, npgfa(iset)
624 rs_grid => rs_rho(igrid_level)
625 map_it(ipgf) =
map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
629 IF (any(map_it(1:npgfa(iset))))
THEN
630 sgfa = first_sgfa(1, iset)
631 ncoa = npgfa(iset)*
ncoset(la_max(iset))
632 m1 = sgfa + nsgfa(iset) - 1
633 ALLOCATE (work(nsgfa(iset), 1))
634 work(1:nsgfa(iset), 1) = aci(sgfa:m1)
637 CALL dgemm(
"N",
"N", ncoa, 1, nsgfa(iset), 1.0_dp, lri_basis_set%sphi(1, sgfa), &
638 SIZE(lri_basis_set%sphi, 1), work(1, 1),
SIZE(work, 1), 0.0_dp, pab(1, 1), &
641 DO ipgf = 1, npgfa(iset)
642 na1 = (ipgf - 1)*
ncoset(la_max(iset))
644 rs_grid => rs_rho(igrid_level)
645 IF (map_it(ipgf))
THEN
647 lb_min=0, lb_max=0, &
648 ra=ra, rb=ra, rp=ra, &
649 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
650 prefactor=1.0_dp, cutoff=1.0_dp)
653 zeta=zeta(ipgf, iset), &
654 la_min=la_min(iset), &
655 lb_max=0, zetb=0.0_dp, lb_min=0, &
656 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
658 pab=pab, o1=na1, o2=0, &
674 IF (exact_1c_terms)
THEN
679 maxsgf_set=maxsgf_set, &
681 ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set))
683 DO ikind = 1,
SIZE(atomic_kind_set)
684 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
685 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=
"ORB")
687 lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
688 sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
690 atom_a = atom_list(iatom)
691 ra(:) =
pbc(particle_set(atom_a)%r, cell)
692 CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, block=p_block, found=found)
693 m1 = maxval(npgfa(1:nseta))
694 ALLOCATE (map_it2(m1, m1))
699 DO ipgf = 1, npgfa(iset)
700 DO jpgf = 1, npgfa(jset)
701 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
703 rs_grid => rs_rho(igrid_level)
704 map_it2(ipgf, jpgf) =
map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
709 IF (any(map_it2(1:npgfa(iset), 1:npgfa(jset))))
THEN
710 ncoa = npgfa(iset)*
ncoset(la_max(iset))
711 sgfa = first_sgfa(1, iset)
712 ncob = npgfa(jset)*
ncoset(la_max(jset))
713 sgfb = first_sgfa(1, jset)
715 CALL dgemm(
"N",
"N", ncoa, nsgfa(jset), nsgfa(iset), &
716 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
717 p_block(sgfa, sgfb),
SIZE(p_block, 1), &
718 0.0_dp, work(1, 1), maxco)
719 CALL dgemm(
"N",
"T", ncoa, ncob, nsgfa(jset), &
720 1.0_dp, work(1, 1), maxco, &
721 sphi_a(1, sgfb),
SIZE(sphi_a, 1), &
722 0.0_dp, pab(1, 1), maxco)
723 DO ipgf = 1, npgfa(iset)
724 DO jpgf = 1, npgfa(jset)
725 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
727 rs_grid => rs_rho(igrid_level)
729 na1 = (ipgf - 1)*
ncoset(la_max(iset))
730 nb1 = (jpgf - 1)*
ncoset(la_max(jset))
732 IF (map_it2(ipgf, jpgf))
THEN
734 la_max=la_max(iset), &
735 lb_min=la_min(jset), &
736 lb_max=la_max(jset), &
737 ra=ra, rb=ra, rp=ra, &
738 zetp=zetp, eps=eps_rho_rspace, &
739 prefactor=1.0_dp, cutoff=1.0_dp)
742 la_max(iset), zeta(ipgf, iset), la_min(iset), &
743 la_max(jset), zeta(jpgf, jset), la_min(jset), &
744 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, na1, nb1, &
757 DEALLOCATE (pab, work)
763 DO igrid_level = 1, gridlevel_info%ngrid_levels
764 CALL pw_zero(mgrid_rspace(igrid_level))
766 pw=mgrid_rspace(igrid_level))
769 DO igrid_level = 1, gridlevel_info%ngrid_levels
770 CALL pw_zero(mgrid_gspace(igrid_level))
772 mgrid_gspace(igrid_level))
773 CALL pw_axpy(mgrid_gspace(igrid_level), lri_rho_g)
782 CALL timestop(handle)
795 SUBROUTINE calculate_rho_core_r3d_rs (rho_core, total_rho, qs_env, calpha, ccore, only_nopaw)
798 REAL(KIND=
dp),
INTENT(OUT) :: total_rho
800 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL :: calpha, ccore
801 LOGICAL,
INTENT(IN),
OPTIONAL :: only_nopaw
803 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_core'
805 INTEGER :: atom_a, handle, iatom, ikind, ithread, &
806 j, natom, npme, nthread, &
808 INTEGER,
DIMENSION(:),
POINTER :: atom_list, cores
809 LOGICAL :: my_only_nopaw, paw_atom
810 REAL(kind=
dp) :: alpha, eps_rho_rspace, radius
811 REAL(kind=
dp),
DIMENSION(3) :: ra
812 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
820 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
823 CALL timeset(routinen, handle)
824 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
825 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
828 my_only_nopaw = .false.
829 IF (
PRESENT(only_nopaw)) my_only_nopaw = only_nopaw
830 IF (
PRESENT(calpha))
THEN
831 cpassert(
PRESENT(ccore))
835 atomic_kind_set=atomic_kind_set, &
836 qs_kind_set=qs_kind_set, &
838 dft_control=dft_control, &
839 particle_set=particle_set, &
841 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
842 auxbas_pw_pool=auxbas_pw_pool)
846 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
848 DO ikind = 1,
SIZE(atomic_kind_set)
849 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
850 IF (
PRESENT(calpha))
THEN
851 alpha = calpha(ikind)
852 pab(1, 1) = ccore(ikind)
854 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
855 IF (my_only_nopaw .AND. paw_atom) cycle
856 CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha, &
857 ccore_charge=pab(1, 1))
860 IF (my_only_nopaw .AND. paw_atom) cycle
861 IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
871 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
873 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
883 IF (npme .GT. 0)
THEN
887 atom_a = atom_list(iatom)
888 ra(:) =
pbc(particle_set(atom_a)%r, cell)
891 lb_min=0, lb_max=0, &
892 ra=ra, rb=ra, rp=ra, &
893 zetp=alpha, eps=eps_rho_rspace, &
894 pab=pab, o1=0, o2=0, &
895 prefactor=-1.0_dp, cutoff=0.0_dp)
898 (/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
900 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
907 IF (
ASSOCIATED(cores))
THEN
912 CALL auxbas_pw_pool%create_pw(rhoc_r)
920 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
922 CALL timestop(handle)
924 END SUBROUTINE calculate_rho_core_r3d_rs
934 SUBROUTINE calculate_rho_core_c1d_gs (rho_core, total_rho, qs_env, calpha, ccore, only_nopaw)
937 REAL(KIND=
dp),
INTENT(OUT) :: total_rho
939 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL :: calpha, ccore
940 LOGICAL,
INTENT(IN),
OPTIONAL :: only_nopaw
942 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_core'
944 INTEGER :: atom_a, handle, iatom, ikind, ithread, &
945 j, natom, npme, nthread, &
947 INTEGER,
DIMENSION(:),
POINTER :: atom_list, cores
948 LOGICAL :: my_only_nopaw, paw_atom
949 REAL(kind=
dp) :: alpha, eps_rho_rspace, radius
950 REAL(kind=
dp),
DIMENSION(3) :: ra
951 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
959 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
962 CALL timeset(routinen, handle)
963 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
964 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
967 my_only_nopaw = .false.
968 IF (
PRESENT(only_nopaw)) my_only_nopaw = only_nopaw
969 IF (
PRESENT(calpha))
THEN
970 cpassert(
PRESENT(ccore))
974 atomic_kind_set=atomic_kind_set, &
975 qs_kind_set=qs_kind_set, &
977 dft_control=dft_control, &
978 particle_set=particle_set, &
980 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
981 auxbas_pw_pool=auxbas_pw_pool)
985 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
987 DO ikind = 1,
SIZE(atomic_kind_set)
988 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
989 IF (
PRESENT(calpha))
THEN
990 alpha = calpha(ikind)
991 pab(1, 1) = ccore(ikind)
993 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
994 IF (my_only_nopaw .AND. paw_atom) cycle
995 CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha, &
996 ccore_charge=pab(1, 1))
999 IF (my_only_nopaw .AND. paw_atom) cycle
1000 IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
1010 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
1012 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
1022 IF (npme .GT. 0)
THEN
1026 atom_a = atom_list(iatom)
1027 ra(:) =
pbc(particle_set(atom_a)%r, cell)
1028 subpatch_pattern = 0
1030 lb_min=0, lb_max=0, &
1031 ra=ra, rb=ra, rp=ra, &
1032 zetp=alpha, eps=eps_rho_rspace, &
1033 pab=pab, o1=0, o2=0, &
1034 prefactor=-1.0_dp, cutoff=0.0_dp)
1037 (/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
1039 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1046 IF (
ASSOCIATED(cores))
THEN
1051 CALL auxbas_pw_pool%create_pw(rhoc_r)
1059 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1061 CALL timestop(handle)
1063 END SUBROUTINE calculate_rho_core_c1d_gs
1078 INTEGER,
INTENT(IN) :: beta, lambda
1080 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_drho_core'
1082 INTEGER :: atom_a, dabqadb_func, handle, iatom, &
1083 ikind, ithread, j, natom, npme, &
1084 nthread, subpatch_pattern
1085 INTEGER,
DIMENSION(:),
POINTER :: atom_list, cores
1086 REAL(kind=
dp) :: alpha, eps_rho_rspace, radius
1087 REAL(kind=
dp),
DIMENSION(3) :: ra
1088 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
1096 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1099 CALL timeset(routinen, handle)
1100 NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
1101 atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
1102 ALLOCATE (pab(1, 1))
1105 atomic_kind_set=atomic_kind_set, &
1106 qs_kind_set=qs_kind_set, &
1108 dft_control=dft_control, &
1109 particle_set=particle_set, &
1111 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1112 auxbas_pw_pool=auxbas_pw_pool)
1116 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1126 cpabort(
"invalid beta")
1128 DO ikind = 1,
SIZE(atomic_kind_set)
1129 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
1131 alpha_core_charge=alpha, ccore_charge=pab(1, 1))
1133 IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
1143 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
1145 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
1155 IF (npme .GT. 0)
THEN
1159 atom_a = atom_list(iatom)
1160 IF (atom_a /= lambda) cycle
1161 ra(:) =
pbc(particle_set(atom_a)%r, cell)
1162 subpatch_pattern = 0
1164 lb_min=0, lb_max=0, &
1165 ra=ra, rb=ra, rp=ra, &
1166 zetp=alpha, eps=eps_rho_rspace, &
1167 pab=pab, o1=0, o2=0, &
1168 prefactor=-1.0_dp, cutoff=0.0_dp)
1171 (/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
1172 radius=radius, ga_gb_function=dabqadb_func, &
1173 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1180 IF (
ASSOCIATED(cores))
THEN
1185 CALL auxbas_pw_pool%create_pw(rhoc_r)
1191 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1193 CALL timestop(handle)
1210 INTEGER,
INTENT(IN) :: iatom_in
1212 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_single_gaussian'
1214 INTEGER :: atom_a, handle, iatom, npme, &
1216 REAL(kind=
dp) :: eps_rho_rspace, radius
1217 REAL(kind=
dp),
DIMENSION(3) :: ra
1218 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
1226 CALL timeset(routinen, handle)
1227 NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool)
1229 ALLOCATE (pab(1, 1))
1233 dft_control=dft_control, &
1235 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1236 auxbas_pw_pool=auxbas_pw_pool)
1239 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1245 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
1246 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
1253 IF (npme .GT. 0)
THEN
1254 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1255 ra(:) =
pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
1256 subpatch_pattern = 0
1258 lb_min=0, lb_max=0, &
1259 ra=ra, rb=ra, rp=ra, &
1260 zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
1261 eps=eps_rho_rspace, &
1262 pab=pab, o1=0, o2=0, &
1263 prefactor=1.0_dp, cutoff=0.0_dp)
1266 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
1268 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1273 CALL auxbas_pw_pool%create_pw(rhoc_r)
1279 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1281 CALL timestop(handle)
1299 REAL(kind=
dp),
DIMENSION(:),
POINTER :: coeff
1300 REAL(kind=
dp),
INTENT(OUT),
OPTIONAL :: total_rho_metal
1303 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_metal'
1305 INTEGER :: atom_a, handle, iatom, j, natom, npme, &
1307 INTEGER,
DIMENSION(:),
POINTER :: cores
1308 REAL(kind=
dp) :: eps_rho_rspace, radius
1309 REAL(kind=
dp),
DIMENSION(3) :: ra
1310 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
1318 CALL timeset(routinen, handle)
1320 NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, cores)
1322 ALLOCATE (pab(1, 1))
1326 dft_control=dft_control, &
1328 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1329 auxbas_pw_pool=auxbas_pw_pool)
1332 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1335 natom =
SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1342 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
1343 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
1353 IF (npme .GT. 0)
THEN
1356 atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1357 ra(:) =
pbc(qs_env%qmmm_env_qm%image_charge_pot%particles_all(atom_a)%r, cell)
1358 subpatch_pattern = 0
1360 lb_min=0, lb_max=0, &
1361 ra=ra, rb=ra, rp=ra, &
1362 zetp=qs_env%qmmm_env_qm%image_charge_pot%eta, &
1363 eps=eps_rho_rspace, &
1364 pab=pab, o1=0, o2=0, &
1365 prefactor=coeff(iatom), cutoff=0.0_dp)
1368 0, qs_env%qmmm_env_qm%image_charge_pot%eta, &
1369 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
1371 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1375 DEALLOCATE (pab, cores)
1377 CALL auxbas_pw_pool%create_pw(rhoc_r)
1381 IF (
PRESENT(total_rho_metal)) &
1386 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1388 CALL timestop(handle)
1406 REAL(kind=
dp),
INTENT(IN) :: eta
1407 INTEGER,
INTENT(IN) :: iatom_in
1409 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_resp_single'
1411 INTEGER :: handle, iatom, npme, subpatch_pattern
1412 REAL(kind=
dp) :: eps_rho_rspace, radius
1413 REAL(kind=
dp),
DIMENSION(3) :: ra
1414 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
1423 CALL timeset(routinen, handle)
1424 NULLIFY (cell, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1427 ALLOCATE (pab(1, 1))
1431 dft_control=dft_control, &
1432 particle_set=particle_set, &
1434 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1435 auxbas_pw_pool=auxbas_pw_pool)
1438 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1444 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
1445 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
1452 IF (npme .GT. 0)
THEN
1453 ra(:) =
pbc(particle_set(iatom)%r, cell)
1454 subpatch_pattern = 0
1456 lb_min=0, lb_max=0, &
1457 ra=ra, rb=ra, rp=ra, &
1458 zetp=eta, eps=eps_rho_rspace, &
1459 pab=pab, o1=0, o2=0, &
1460 prefactor=1.0_dp, cutoff=0.0_dp)
1463 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, rs_rho, &
1465 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1470 CALL auxbas_pw_pool%create_pw(rhoc_r)
1476 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1478 CALL timestop(handle)
1494 SUBROUTINE calculate_rho_resp_all_r3d_rs (rho_resp, coeff, natom, eta, qs_env)
1497 REAL(KIND=
dp),
DIMENSION(:),
POINTER :: coeff
1498 INTEGER,
INTENT(IN) :: natom
1499 REAL(KIND=
dp),
INTENT(IN) :: eta
1502 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_resp_all'
1504 INTEGER :: handle, iatom, j, npme, subpatch_pattern
1505 INTEGER,
DIMENSION(:),
POINTER :: cores
1506 REAL(kind=
dp) :: eps_rho_rspace, radius
1507 REAL(kind=
dp),
DIMENSION(3) :: ra
1508 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
1517 CALL timeset(routinen, handle)
1519 NULLIFY (cell, cores, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1522 ALLOCATE (pab(1, 1))
1526 dft_control=dft_control, &
1527 particle_set=particle_set, &
1529 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1530 auxbas_pw_pool=auxbas_pw_pool)
1533 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1541 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
1542 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
1552 IF (npme .GT. 0)
THEN
1555 ra(:) =
pbc(particle_set(iatom)%r, cell)
1556 subpatch_pattern = 0
1558 lb_min=0, lb_max=0, &
1559 ra=ra, rb=ra, rp=ra, &
1560 zetp=eta, eps=eps_rho_rspace, &
1561 pab=pab, o1=0, o2=0, &
1562 prefactor=coeff(iatom), cutoff=0.0_dp)
1566 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
1568 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1572 DEALLOCATE (pab, cores)
1574 CALL auxbas_pw_pool%create_pw(rhoc_r)
1579 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1581 CALL timestop(handle)
1583 END SUBROUTINE calculate_rho_resp_all_r3d_rs
1596 SUBROUTINE calculate_rho_resp_all_c1d_gs (rho_resp, coeff, natom, eta, qs_env)
1599 REAL(KIND=
dp),
DIMENSION(:),
POINTER :: coeff
1600 INTEGER,
INTENT(IN) :: natom
1601 REAL(KIND=
dp),
INTENT(IN) :: eta
1604 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_resp_all'
1606 INTEGER :: handle, iatom, j, npme, subpatch_pattern
1607 INTEGER,
DIMENSION(:),
POINTER :: cores
1608 REAL(kind=
dp) :: eps_rho_rspace, radius
1609 REAL(kind=
dp),
DIMENSION(3) :: ra
1610 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab
1619 CALL timeset(routinen, handle)
1621 NULLIFY (cell, cores, dft_control, pab, pw_env, rs_rho, auxbas_pw_pool, &
1624 ALLOCATE (pab(1, 1))
1628 dft_control=dft_control, &
1629 particle_set=particle_set, &
1631 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
1632 auxbas_pw_pool=auxbas_pw_pool)
1635 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1643 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed)
THEN
1644 IF (
modulo(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos)
THEN
1654 IF (npme .GT. 0)
THEN
1657 ra(:) =
pbc(particle_set(iatom)%r, cell)
1658 subpatch_pattern = 0
1660 lb_min=0, lb_max=0, &
1661 ra=ra, rb=ra, rp=ra, &
1662 zetp=eta, eps=eps_rho_rspace, &
1663 pab=pab, o1=0, o2=0, &
1664 prefactor=coeff(iatom), cutoff=0.0_dp)
1668 0, 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), coeff(iatom), pab, 0, 0, rs_rho, &
1670 use_subpatch=.true., subpatch_pattern=subpatch_pattern)
1674 DEALLOCATE (pab, cores)
1676 CALL auxbas_pw_pool%create_pw(rhoc_r)
1681 CALL auxbas_pw_pool%give_back_pw(rhoc_r)
1683 CALL timestop(handle)
1685 END SUBROUTINE calculate_rho_resp_all_c1d_gs
1713 ks_env, soft_valid, compute_tau, compute_grad, &
1714 basis_type, der_type, idir, task_list_external, pw_env_external)
1716 TYPE(
dbcsr_type),
OPTIONAL,
TARGET :: matrix_p
1718 POINTER :: matrix_p_kp
1721 REAL(kind=
dp),
INTENT(OUT),
OPTIONAL :: total_rho
1723 LOGICAL,
INTENT(IN),
OPTIONAL :: soft_valid, compute_tau, compute_grad
1724 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: basis_type
1725 INTEGER,
INTENT(IN),
OPTIONAL :: der_type, idir
1727 TYPE(
pw_env_type),
OPTIONAL,
POINTER :: pw_env_external
1729 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_rho_elec'
1731 CHARACTER(LEN=default_string_length) :: my_basis_type
1732 INTEGER :: ga_gb_function, handle, ilevel, img, &
1734 LOGICAL :: any_distributed, my_compute_grad, &
1735 my_compute_tau, my_soft_valid
1736 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_images
1743 CALL timeset(routinen, handle)
1745 NULLIFY (matrix_images, dft_control, pw_env, rs_rho, task_list)
1748 my_compute_tau = .false.
1749 IF (
PRESENT(compute_tau)) my_compute_tau = compute_tau
1750 my_compute_grad = .false.
1751 IF (
PRESENT(compute_grad)) my_compute_grad = compute_grad
1752 IF (
PRESENT(der_type))
THEN
1753 SELECT CASE (der_type)
1775 cpabort(
"Unknown der_type")
1777 ELSE IF (my_compute_tau)
THEN
1779 ELSE IF (my_compute_grad)
THEN
1780 cpassert(
PRESENT(idir))
1789 cpabort(
"invalid idir")
1796 my_basis_type =
"ORB"
1797 IF (
PRESENT(basis_type)) my_basis_type = basis_type
1798 cpassert(my_basis_type ==
"ORB" .OR.
PRESENT(task_list_external))
1801 my_soft_valid = .false.
1802 IF (
PRESENT(soft_valid)) my_soft_valid = soft_valid
1803 IF (
PRESENT(task_list_external))
THEN
1804 task_list => task_list_external
1805 ELSEIF (my_soft_valid)
THEN
1806 CALL get_ks_env(ks_env, task_list_soft=task_list)
1810 cpassert(
ASSOCIATED(task_list))
1813 IF (
PRESENT(pw_env_external))
THEN
1814 pw_env => pw_env_external
1818 cpassert(
ASSOCIATED(pw_env))
1822 nlevels =
SIZE(rs_rho)
1823 group = rs_rho(1)%desc%group
1826 any_distributed = .false.
1827 DO ilevel = 1, nlevels
1828 any_distributed = any_distributed .OR. rs_rho(ilevel)%desc%distributed
1832 CALL get_ks_env(ks_env, dft_control=dft_control)
1833 nimages = dft_control%nimages
1834 ALLOCATE (matrix_images(nimages))
1835 IF (
PRESENT(matrix_p_kp))
THEN
1836 cpassert(.NOT.
PRESENT(matrix_p))
1838 matrix_images(img)%matrix => matrix_p_kp(img)%matrix
1841 cpassert(
PRESENT(matrix_p) .AND. nimages == 1)
1842 matrix_images(1)%matrix => matrix_p
1846 IF (any_distributed)
THEN
1851 DEALLOCATE (matrix_images)
1855 ga_gb_function=ga_gb_function, &
1856 pab_blocks=task_list%pab_buffer, &
1863 CALL timestop(handle)
1880 soft_valid, basis_type)
1882 TYPE(
dbcsr_type),
OPTIONAL,
TARGET :: matrix_p
1884 POINTER :: matrix_p_kp
1888 LOGICAL,
INTENT(IN),
OPTIONAL :: soft_valid
1889 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: basis_type
1891 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_drho_elec'
1893 CHARACTER(LEN=default_string_length) :: my_basis_type
1894 INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, idir, igrid_level, ikind, &
1895 ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
1896 jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
1897 ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
1898 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1900 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb
1901 LOGICAL :: atom_pair_changed, distributed_rs_grids, &
1902 do_kp, found, my_soft, use_subpatch
1903 REAL(kind=
dp) :: eps_rho_rspace, f, prefactor, radius, &
1905 REAL(kind=
dp),
DIMENSION(3) :: ra, rab, rab_inv, rb, rp
1906 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: p_block, pab, sphi_a, sphi_b, work, &
1908 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: pabt, workt
1909 TYPE(
atom_pair_type),
DIMENSION(:),
POINTER :: atom_pair_recv, atom_pair_send
1919 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1924 TYPE(
task_type),
DIMENSION(:),
POINTER :: tasks
1926 CALL timeset(routinen, handle)
1928 cpassert(
PRESENT(matrix_p) .OR.
PRESENT(matrix_p_kp))
1929 do_kp =
PRESENT(matrix_p_kp)
1931 NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
1932 sab_orb, particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, &
1933 lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, &
1934 sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
1938 IF (
PRESENT(soft_valid)) my_soft = soft_valid
1940 IF (
PRESENT(basis_type))
THEN
1941 my_basis_type = basis_type
1943 my_basis_type =
"ORB"
1947 qs_kind_set=qs_kind_set, &
1949 dft_control=dft_control, &
1950 particle_set=particle_set, &
1954 SELECT CASE (my_basis_type)
1957 task_list=task_list, &
1958 task_list_soft=task_list_soft)
1961 task_list_soft=task_list_soft)
1962 CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
1966 gridlevel_info => pw_env%gridlevel_info
1972 maxsgf_set=maxsgf_set, &
1973 basis_type=my_basis_type)
1974 CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
1975 CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
1978 nimages = dft_control%nimages
1979 cpassert(nimages == 1 .OR. do_kp)
1981 natoms =
SIZE(particle_set)
1984 IF (my_soft) task_list => task_list_soft
1985 cpassert(
ASSOCIATED(task_list))
1986 tasks => task_list%tasks
1987 atom_pair_send => task_list%atom_pair_send
1988 atom_pair_recv => task_list%atom_pair_recv
1989 ntasks = task_list%ntasks
1992 cpassert(
ASSOCIATED(pw_env))
1993 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
1994 DO igrid_level = 1, gridlevel_info%ngrid_levels
1995 distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
1998 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2003 ALLOCATE (deltap(nimages))
2004 IF (distributed_rs_grids)
THEN
2011 CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
2015 CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name=
"DeltaP")
2020 deltap(img)%matrix => matrix_p_kp(img)%matrix
2023 deltap(1)%matrix => matrix_p
2028 IF (distributed_rs_grids)
THEN
2030 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
2031 nimages=nimages, scatter=.true.)
2037 pab => pabt(:, :, ithread)
2038 work => workt(:, :, ithread)
2040 loop_xyz:
DO idir = 1, 3
2042 DO igrid_level = 1, gridlevel_info%ngrid_levels
2046 iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
2047 ikind_old = -1; jkind_old = -1; img_old = -1
2048 loop_tasks:
DO itask = 1, ntasks
2051 igrid_level = tasks(itask)%grid_level
2052 img = tasks(itask)%image
2053 iatom = tasks(itask)%iatom
2054 jatom = tasks(itask)%jatom
2055 iset = tasks(itask)%iset
2056 jset = tasks(itask)%jset
2057 ipgf = tasks(itask)%ipgf
2058 jpgf = tasks(itask)%jpgf
2060 ikind = particle_set(iatom)%atomic_kind%kind_number
2061 jkind = particle_set(jatom)%atomic_kind%kind_number
2063 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old .OR. img .NE. img_old)
THEN
2065 IF (iatom .NE. iatom_old) ra(:) =
pbc(particle_set(iatom)%r, cell)
2067 IF (iatom <= jatom)
THEN
2075 IF (ikind .NE. ikind_old)
THEN
2077 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2078 basis_type=
"ORB_SOFT")
2080 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2081 basis_type=my_basis_type)
2084 first_sgf=first_sgfa, &
2094 IF (jkind .NE. jkind_old)
THEN
2096 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2097 basis_type=
"ORB_SOFT")
2099 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2100 basis_type=my_basis_type)
2103 first_sgf=first_sgfb, &
2114 row=brow, col=bcol, block=p_block, found=found)
2122 atom_pair_changed = .true.
2126 atom_pair_changed = .false.
2130 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset)
THEN
2132 ncoa = npgfa(iset)*
ncoset(la_max(iset))
2133 sgfa = first_sgfa(1, iset)
2134 ncob = npgfb(jset)*
ncoset(lb_max(jset))
2135 sgfb = first_sgfb(1, jset)
2137 IF (iatom <= jatom)
THEN
2138 CALL dgemm(
"N",
"N", ncoa, nsgfb(jset), nsgfa(iset), &
2139 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
2140 p_block(sgfa, sgfb),
SIZE(p_block, 1), &
2141 0.0_dp, work(1, 1), maxco)
2142 CALL dgemm(
"N",
"T", ncoa, ncob, nsgfb(jset), &
2143 1.0_dp, work(1, 1), maxco, &
2144 sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
2145 0.0_dp, pab(1, 1), maxco)
2147 CALL dgemm(
"N",
"N", ncob, nsgfa(iset), nsgfb(jset), &
2148 1.0_dp, sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
2149 p_block(sgfb, sgfa),
SIZE(p_block, 1), &
2150 0.0_dp, work(1, 1), maxco)
2151 CALL dgemm(
"N",
"T", ncob, ncoa, nsgfa(iset), &
2152 1.0_dp, work(1, 1), maxco, &
2153 sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
2154 0.0_dp, pab(1, 1), maxco)
2162 rab(:) = tasks(itask)%rab
2163 rb(:) = ra(:) + rab(:)
2164 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
2166 f = zetb(jpgf, jset)/zetp
2167 rp(:) = ra(:) + f*rab(:)
2168 prefactor = exp(-zeta(ipgf, iset)*f*dot_product(rab, rab))
2170 lb_min=lb_min(jset), lb_max=lb_max(jset), &
2171 ra=ra, rb=rb, rp=rp, &
2172 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
2173 prefactor=prefactor, cutoff=1.0_dp)
2175 na1 = (ipgf - 1)*
ncoset(la_max(iset)) + 1
2176 na2 = ipgf*
ncoset(la_max(iset))
2177 nb1 = (jpgf - 1)*
ncoset(lb_max(jset)) + 1
2178 nb2 = jpgf*
ncoset(lb_max(jset))
2181 IF (iatom == jatom .AND. img == 1)
THEN
2188 IF (rs_rho(igrid_level)%desc%distributed)
THEN
2190 IF (tasks(itask)%dist_type .EQ. 2)
THEN
2191 use_subpatch = .true.
2193 use_subpatch = .false.
2196 use_subpatch = .false.
2207 cpabort(
"invalid idir")
2210 IF (iatom <= jatom)
THEN
2212 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2213 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2214 ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2215 rs_rho(igrid_level), &
2216 radius=radius, ga_gb_function=dabqadb_func, &
2217 use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
2221 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2222 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2223 rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2224 rs_rho(igrid_level), &
2225 radius=radius, ga_gb_function=dabqadb_func, &
2226 use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
2231 CALL density_rs2pw(pw_env, rs_rho, drho(idir), drho_gspace(idir))
2236 IF (distributed_rs_grids)
THEN
2240 NULLIFY (deltap(img)%matrix)
2245 DEALLOCATE (pabt, workt)
2247 CALL timestop(handle)
2267 soft_valid, basis_type, beta, lambda)
2269 TYPE(
dbcsr_type),
OPTIONAL,
TARGET :: matrix_p
2271 POINTER :: matrix_p_kp
2275 LOGICAL,
INTENT(IN),
OPTIONAL :: soft_valid
2276 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: basis_type
2277 INTEGER,
INTENT(IN) :: beta, lambda
2279 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_drho_elec_dR'
2281 CHARACTER(LEN=default_string_length) :: my_basis_type
2282 INTEGER :: bcol, brow, dabqadb_func, handle, iatom, iatom_old, igrid_level, ikind, &
2283 ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
2284 jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
2285 ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
2286 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
2288 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb
2289 LOGICAL :: atom_pair_changed, distributed_rs_grids, &
2290 do_kp, found, my_soft, use_subpatch
2291 REAL(kind=
dp) :: eps_rho_rspace, f, prefactor, radius, &
2293 REAL(kind=
dp),
DIMENSION(3) :: ra, rab, rab_inv, rb, rp
2294 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: p_block, pab, sphi_a, sphi_b, work, &
2296 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: pabt, workt
2297 TYPE(
atom_pair_type),
DIMENSION(:),
POINTER :: atom_pair_recv, atom_pair_send
2305 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2310 TYPE(
task_type),
DIMENSION(:),
POINTER :: tasks
2312 CALL timeset(routinen, handle)
2314 cpassert(
PRESENT(matrix_p) .OR.
PRESENT(matrix_p_kp))
2315 do_kp =
PRESENT(matrix_p_kp)
2317 NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
2318 particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, lb_max, &
2319 lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, sphi_b, &
2320 zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
2324 IF (
PRESENT(soft_valid)) my_soft = soft_valid
2326 IF (
PRESENT(basis_type))
THEN
2327 my_basis_type = basis_type
2329 my_basis_type =
"ORB"
2333 qs_kind_set=qs_kind_set, &
2335 dft_control=dft_control, &
2336 particle_set=particle_set, &
2339 SELECT CASE (my_basis_type)
2342 task_list=task_list, &
2343 task_list_soft=task_list_soft)
2346 task_list_soft=task_list_soft)
2347 CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list)
2351 gridlevel_info => pw_env%gridlevel_info
2357 maxsgf_set=maxsgf_set, &
2358 basis_type=my_basis_type)
2359 CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
2360 CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
2363 nimages = dft_control%nimages
2364 cpassert(nimages == 1 .OR. do_kp)
2366 natoms =
SIZE(particle_set)
2369 IF (my_soft) task_list => task_list_soft
2370 cpassert(
ASSOCIATED(task_list))
2371 tasks => task_list%tasks
2372 atom_pair_send => task_list%atom_pair_send
2373 atom_pair_recv => task_list%atom_pair_recv
2374 ntasks = task_list%ntasks
2377 cpassert(
ASSOCIATED(pw_env))
2378 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
2379 DO igrid_level = 1, gridlevel_info%ngrid_levels
2380 distributed_rs_grids = rs_rho(igrid_level)%desc%distributed
2383 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2388 ALLOCATE (deltap(nimages))
2389 IF (distributed_rs_grids)
THEN
2396 CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
2400 CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name=
"DeltaP")
2405 deltap(img)%matrix => matrix_p_kp(img)%matrix
2408 deltap(1)%matrix => matrix_p
2413 IF (distributed_rs_grids)
THEN
2415 atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
2416 nimages=nimages, scatter=.true.)
2422 pab => pabt(:, :, ithread)
2423 work => workt(:, :, ithread)
2425 DO igrid_level = 1, gridlevel_info%ngrid_levels
2429 iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
2430 ikind_old = -1; jkind_old = -1; img_old = -1
2431 loop_tasks:
DO itask = 1, ntasks
2434 igrid_level = tasks(itask)%grid_level
2435 img = tasks(itask)%image
2436 iatom = tasks(itask)%iatom
2437 jatom = tasks(itask)%jatom
2438 iset = tasks(itask)%iset
2439 jset = tasks(itask)%jset
2440 ipgf = tasks(itask)%ipgf
2441 jpgf = tasks(itask)%jpgf
2443 ikind = particle_set(iatom)%atomic_kind%kind_number
2444 jkind = particle_set(jatom)%atomic_kind%kind_number
2446 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old .OR. img .NE. img_old)
THEN
2448 IF (iatom .NE. iatom_old) ra(:) =
pbc(particle_set(iatom)%r, cell)
2450 IF (iatom <= jatom)
THEN
2458 IF (ikind .NE. ikind_old)
THEN
2460 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2461 basis_type=
"ORB_SOFT")
2463 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
2464 basis_type=my_basis_type)
2467 first_sgf=first_sgfa, &
2477 IF (jkind .NE. jkind_old)
THEN
2479 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2480 basis_type=
"ORB_SOFT")
2482 CALL get_qs_kind(qs_kind_set(jkind), basis_set=orb_basis_set, &
2483 basis_type=my_basis_type)
2486 first_sgf=first_sgfb, &
2497 row=brow, col=bcol, block=p_block, found=found)
2505 atom_pair_changed = .true.
2509 atom_pair_changed = .false.
2513 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset)
THEN
2515 ncoa = npgfa(iset)*
ncoset(la_max(iset))
2516 sgfa = first_sgfa(1, iset)
2517 ncob = npgfb(jset)*
ncoset(lb_max(jset))
2518 sgfb = first_sgfb(1, jset)
2520 IF (iatom <= jatom)
THEN
2521 CALL dgemm(
"N",
"N", ncoa, nsgfb(jset), nsgfa(iset), &
2522 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
2523 p_block(sgfa, sgfb),
SIZE(p_block, 1), &
2524 0.0_dp, work(1, 1), maxco)
2525 CALL dgemm(
"N",
"T", ncoa, ncob, nsgfb(jset), &
2526 1.0_dp, work(1, 1), maxco, &
2527 sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
2528 0.0_dp, pab(1, 1), maxco)
2530 CALL dgemm(
"N",
"N", ncob, nsgfa(iset), nsgfb(jset), &
2531 1.0_dp, sphi_b(1, sgfb),
SIZE(sphi_b, 1), &
2532 p_block(sgfb, sgfa),
SIZE(p_block, 1), &
2533 0.0_dp, work(1, 1), maxco)
2534 CALL dgemm(
"N",
"T", ncob, ncoa, nsgfa(iset), &
2535 1.0_dp, work(1, 1), maxco, &
2536 sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
2537 0.0_dp, pab(1, 1), maxco)
2545 rab(:) = tasks(itask)%rab
2546 rb(:) = ra(:) + rab(:)
2547 zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
2549 f = zetb(jpgf, jset)/zetp
2550 rp(:) = ra(:) + f*rab(:)
2551 prefactor = exp(-zeta(ipgf, iset)*f*dot_product(rab, rab))
2553 lb_min=lb_min(jset), lb_max=lb_max(jset), &
2554 ra=ra, rb=rb, rp=rp, &
2555 zetp=zetp, eps=eps_rho_rspace, &
2556 prefactor=prefactor, cutoff=1.0_dp)
2558 na1 = (ipgf - 1)*
ncoset(la_max(iset)) + 1
2559 na2 = ipgf*
ncoset(la_max(iset))
2560 nb1 = (jpgf - 1)*
ncoset(lb_max(jset)) + 1
2561 nb2 = jpgf*
ncoset(lb_max(jset))
2564 IF (iatom == jatom .AND. img == 1)
THEN
2571 IF (rs_rho(igrid_level)%desc%distributed)
THEN
2573 IF (tasks(itask)%dist_type .EQ. 2)
THEN
2574 use_subpatch = .true.
2576 use_subpatch = .false.
2579 use_subpatch = .false.
2590 cpabort(
"invalid beta")
2593 IF (iatom <= jatom)
THEN
2594 IF (iatom == lambda) &
2596 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2597 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2598 ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2599 rsgrid=rs_rho(igrid_level), &
2600 ga_gb_function=dabqadb_func, radius=radius, &
2601 use_subpatch=use_subpatch, &
2602 subpatch_pattern=tasks(itask)%subpatch_pattern)
2603 IF (jatom == lambda) &
2605 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2606 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2607 ra, rab, scale, pab, na1 - 1, nb1 - 1, &
2608 rsgrid=rs_rho(igrid_level), &
2609 ga_gb_function=dabqadb_func + 3, radius=radius, &
2610 use_subpatch=use_subpatch, &
2611 subpatch_pattern=tasks(itask)%subpatch_pattern)
2614 IF (jatom == lambda) &
2616 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2617 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2618 rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2619 rs_rho(igrid_level), &
2620 ga_gb_function=dabqadb_func, radius=radius, &
2621 use_subpatch=use_subpatch, &
2622 subpatch_pattern=tasks(itask)%subpatch_pattern)
2623 IF (iatom == lambda) &
2625 lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
2626 la_max(iset), zeta(ipgf, iset), la_min(iset), &
2627 rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
2628 rs_rho(igrid_level), &
2629 ga_gb_function=dabqadb_func + 3, radius=radius, &
2630 use_subpatch=use_subpatch, &
2631 subpatch_pattern=tasks(itask)%subpatch_pattern)
2639 IF (distributed_rs_grids)
THEN
2643 NULLIFY (deltap(img)%matrix)
2648 DEALLOCATE (pabt, workt)
2650 CALL timestop(handle)
2673 atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
2674 pw_env, required_function, basis_type)
2679 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2684 INTEGER,
INTENT(IN) :: required_function
2685 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: basis_type
2687 CHARACTER(LEN=*),
PARAMETER :: routinen =
'collocate_single_gaussian'
2689 CHARACTER(LEN=default_string_length) :: my_basis_type
2690 INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
2691 my_index, my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
2692 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: where_is_the_point
2693 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, npgfa, nsgfa
2694 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa
2696 REAL(kind=
dp) :: dab, eps_rho_rspace, radius, scale
2697 REAL(kind=
dp),
DIMENSION(3) :: ra
2698 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab, sphi_a, zeta
2707 IF (
PRESENT(basis_type))
THEN
2708 my_basis_type = basis_type
2710 my_basis_type =
"ORB"
2713 CALL timeset(routinen, handle)
2715 NULLIFY (orb_basis_set, pab, la_max, la_min, npgfa, nsgfa, sphi_a, &
2716 zeta, first_sgfa, rs_rho, pw_pools)
2719 cpassert(
ASSOCIATED(pw_env))
2720 CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
2721 gridlevel_info=gridlevel_info)
2727 DO igrid_level = 1, gridlevel_info%ngrid_levels
2731 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2736 maxsgf_set=maxsgf_set, &
2737 basis_type=my_basis_type)
2739 ALLOCATE (pab(maxco, 1))
2742 group = mgrid_rspace(1)%pw_grid%para%group
2743 my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
2744 group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
2745 ALLOCATE (where_is_the_point(0:group_size - 1))
2748 ikind = particle_set(iatom)%atomic_kind%kind_number
2749 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
2751 first_sgf=first_sgfa, &
2759 ra(:) =
pbc(particle_set(iatom)%r, cell)
2764 ncoa = npgfa(iset)*
ncoset(la_max(iset))
2765 sgfa = first_sgfa(1, iset)
2769 DO i = 1, nsgfa(iset)
2770 IF (offset + i == required_function)
THEN
2779 pab(1:ncoa, 1) = sphi_a(1:ncoa, sgfa + my_index - 1)
2781 DO ipgf = 1, npgfa(iset)
2783 na1 = (ipgf - 1)*
ncoset(la_max(iset)) + 1
2784 na2 = ipgf*
ncoset(la_max(iset))
2789 IF (
map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos))
THEN
2791 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2792 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
2793 prefactor=1.0_dp, cutoff=1.0_dp)
2797 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
2798 scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
2806 offset = offset + nsgfa(iset)
2812 DO igrid_level = 1, gridlevel_info%ngrid_levels
2814 mgrid_rspace(igrid_level))
2818 DO igrid_level = 1, gridlevel_info%ngrid_levels
2820 mgrid_gspace(igrid_level))
2821 CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
2833 CALL timestop(handle)
2856 atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, &
2859 INTEGER,
INTENT(IN) :: ivector
2863 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2868 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: basis_type
2870 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_wavefunction'
2872 INTEGER :: handle, i, nao
2874 REAL(kind=
dp) :: eps_rho_rspace
2875 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigenvector
2877 CALL timeset(routinen, handle)
2880 ALLOCATE (eigenvector(nao))
2885 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
2888 atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
2889 eps_rho_rspace, basis_type)
2891 DEALLOCATE (eigenvector)
2893 CALL timestop(handle)
2918 atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
2919 eps_rho_rspace, basis_type)
2920 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: vector
2924 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2928 REAL(kind=
dp),
INTENT(IN) :: eps_rho_rspace
2929 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: basis_type
2931 CHARACTER(LEN=*),
PARAMETER :: routinen =
'collocate_function'
2933 CHARACTER(LEN=default_string_length) :: my_basis_type
2934 INTEGER :: group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, maxsgf_set, &
2935 my_pos, na1, na2, natom, ncoa, nseta, offset, sgfa
2936 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: where_is_the_point
2937 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, npgfa, nsgfa
2938 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa
2939 REAL(kind=
dp) :: dab, radius, scale
2940 REAL(kind=
dp),
DIMENSION(3) :: ra
2941 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pab, sphi_a, work, zeta
2950 CALL timeset(routinen, handle)
2952 IF (
PRESENT(basis_type))
THEN
2953 my_basis_type = basis_type
2955 my_basis_type =
"ORB"
2958 NULLIFY (orb_basis_set, pab, work, la_max, la_min, &
2959 npgfa, nsgfa, sphi_a, zeta, first_sgfa, rs_rho, pw_pools)
2962 cpassert(
ASSOCIATED(pw_env))
2963 CALL pw_env_get(pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
2964 gridlevel_info=gridlevel_info)
2970 DO igrid_level = 1, gridlevel_info%ngrid_levels
2978 maxsgf_set=maxsgf_set, &
2979 basis_type=my_basis_type)
2981 ALLOCATE (pab(maxco, 1))
2982 ALLOCATE (work(maxco, 1))
2985 group = mgrid_rspace(1)%pw_grid%para%group
2986 my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
2987 group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
2988 ALLOCATE (where_is_the_point(0:group_size - 1))
2991 ikind = particle_set(iatom)%atomic_kind%kind_number
2992 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=my_basis_type)
2994 first_sgf=first_sgfa, &
3002 ra(:) =
pbc(particle_set(iatom)%r, cell)
3007 ncoa = npgfa(iset)*
ncoset(la_max(iset))
3008 sgfa = first_sgfa(1, iset)
3010 DO i = 1, nsgfa(iset)
3011 work(i, 1) = vector(offset + i)
3014 CALL dgemm(
"N",
"N", ncoa, 1, nsgfa(iset), &
3015 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
3016 work(1, 1),
SIZE(work, 1), &
3017 0.0_dp, pab(1, 1),
SIZE(pab, 1))
3019 DO ipgf = 1, npgfa(iset)
3021 na1 = (ipgf - 1)*
ncoset(la_max(iset)) + 1
3022 na2 = ipgf*
ncoset(la_max(iset))
3027 IF (
map_gaussian_here(rs_rho(igrid_level), cell%h_inv, ra, offset, group_size, my_pos))
THEN
3029 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
3030 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
3031 prefactor=1.0_dp, cutoff=1.0_dp)
3035 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
3036 scale, pab, na1 - 1, 0, rs_rho(igrid_level), &
3042 offset = offset + nsgfa(iset)
3048 DO igrid_level = 1, gridlevel_info%ngrid_levels
3050 mgrid_rspace(igrid_level))
3054 DO igrid_level = 1, gridlevel_info%ngrid_levels
3056 mgrid_gspace(igrid_level))
3057 CALL pw_axpy(mgrid_gspace(igrid_level), rho_gspace)
3070 CALL timestop(handle)
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
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.
All kind of helpful little routines.
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
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_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
DBCSR operations in CP2K.
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_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition of the atomic potential types.
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Fortran API for the grid package, which is written in C.
integer, parameter, public grid_func_core_x
integer, parameter, public grid_func_dab_z
subroutine, public grid_collocate_task_list(task_list, ga_gb_function, pab_blocks, rs_grids)
Collocate all tasks of in given list onto given grids.
integer, parameter, public grid_func_dzdx
integer, parameter, public grid_func_dzdz
integer, parameter, public grid_func_dydz
integer, parameter, public grid_func_dxdy
integer, parameter, public grid_func_dabpadb_y
integer, parameter, public grid_func_dab_y
integer, parameter, public grid_func_dxdx
integer, parameter, public grid_func_dadb
integer, parameter, public grid_func_dydy
integer, parameter, public grid_func_dabpadb_z
integer, parameter, public grid_func_dabpadb_x
integer, parameter, public grid_func_dx
integer, parameter, public grid_func_dz
integer, parameter, public grid_func_ab
integer, parameter, public grid_func_core_y
integer, parameter, public grid_func_dab_x
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
integer, parameter, public grid_func_core_z
integer, parameter, public grid_func_dy
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Utility routines for the memory handling.
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :, :), allocatable, public coset
Define the data structure for the particle information.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_drho_elec_dr(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, soft_valid, basis_type, beta, lambda)
Computes the gradient wrt. nuclear coordinates of a density on the grid The density is given in terms...
subroutine, public calculate_drho_elec(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, soft_valid, basis_type)
computes the gradient of the density corresponding to a given density matrix on the grid
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 collocate_function(vector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, eps_rho_rspace, basis_type)
maps a given function on the grid
subroutine, public calculate_rho_metal(rho_metal, coeff, total_rho_metal, qs_env)
computes the image charge density on the grid (including coeffcients)
subroutine, public calculate_rho_resp_single(rho_gb, qs_env, eta, iatom_in)
collocate a single Gaussian on the grid for periodic RESP fitting
subroutine, public calculate_rho_nlcc(rho_nlcc, qs_env)
computes the density of the non-linear core correction on the grid
subroutine, public calculate_lri_rho_elec(lri_rho_g, lri_rho_r, qs_env, lri_coef, total_rho, basis_type, exact_1c_terms, pmat, atomlist)
Collocates the fitted lri density on a grid.
subroutine, public collocate_single_gaussian(rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, required_function, basis_type)
maps a single gaussian on the grid
subroutine, public calculate_rho_single_gaussian(rho_gb, qs_env, iatom_in)
collocate a single Gaussian on the grid
subroutine, public calculate_ppl_grid(vppl, qs_env)
computes the local pseudopotential (without erf term) on the grid
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
subroutine, public calculate_drho_core(drho_core, qs_env, beta, lambda)
Computes the derivative of the density of the core charges with respect to the nuclear coordinates on...
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.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
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)
...
Define the neighbor list data types and the corresponding functionality.
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
subroutine, public transfer_rs2pw(rs, pw)
...
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Transfers densities from PW to RS grids and potentials from PW to RS.
subroutine, public density_rs2pw(pw_env, rs_rho, rho, rho_gspace)
given partial densities on the realspace multigrids, computes the full density on the plane wave grid...
generate the tasks lists used by collocate and integrate routines
subroutine, public rs_scatter_matrices(src_matrices, dest_buffer, task_list, group)
Scatters dbcsr matrix blocks and receives them into a buffer as needed before collocation.
subroutine, public rs_distribute_matrix(rs_descs, pmats, atom_pair_send, atom_pair_recv, nimages, scatter, hmats)
redistributes the matrix so that it can be used in realspace operations i.e. according to the task li...
subroutine, public rs_copy_to_buffer(src_matrices, dest_buffer, task_list)
Copies the DBCSR blocks into buffer, replaces rs_scatter_matrix for non-distributed grids.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
contained for different pw related things
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...