97 SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
98 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
99 nimages, cell_to_index, basis_type, deltaR, atcore)
101 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_h, matrix_p
104 LOGICAL,
INTENT(IN) :: calculate_forces
105 LOGICAL :: use_virial
107 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
111 POINTER :: sab_orb, sac_ppl
112 INTEGER,
INTENT(IN) :: nimages
113 INTEGER,
DIMENSION(:, :, :),
OPTIONAL,
POINTER :: cell_to_index
114 CHARACTER(LEN=*),
INTENT(IN) :: basis_type
115 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
117 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT), &
120 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_core_ppl'
121 INTEGER,
PARAMETER :: nexp_max = 30
123 INTEGER :: atom_a, handle, i, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, &
124 katom, kkind, ldai, ldsab, maxco, maxder, maxl, maxlgto, maxlppl, maxnset, maxsgf, mepos, &
125 n_local, natom, ncoa, ncob, nexp_lpot, nexp_ppl, nkind, nloc, nseta, nsetb, nthread, &
126 sgfa, sgfb, slmax, slot
127 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
128 INTEGER,
DIMENSION(0:10) :: npot
129 INTEGER,
DIMENSION(1:10) :: nrloc
130 INTEGER,
DIMENSION(1:15, 0:10) :: nrpot
131 INTEGER,
DIMENSION(3) :: cellind
132 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, &
133 nct_lpot, npgfa, npgfb, nsgfa, nsgfb
134 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb
135 INTEGER,
DIMENSION(nexp_max) :: nct_ppl
136 LOGICAL :: do_dr, do_upf, doat, dokp, ecp_local, &
137 ecp_semi_local, found, have_upf, &
138 libgrpp_local, lpotextended, &
140 REAL(kind=
dp) :: alpha, atk0, atk1, dab, dac, dbc, f0, &
142 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: upf_v, work
143 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work, upf_da, &
145 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: hab, pab
146 REAL(kind=
dp),
ALLOCATABLE, &
147 DIMENSION(:, :, :, :, :) :: hab2
148 REAL(kind=
dp),
DIMENSION(1:10) :: aloc, bloc
149 REAL(kind=
dp),
DIMENSION(1:15, 0:10) :: apot, bpot
150 REAL(kind=
dp),
DIMENSION(3) :: force_a, force_b, rab, rac, rbc
151 REAL(kind=
dp),
DIMENSION(3, 3) :: pv_thread
153 DIMENSION(:),
POINTER :: ap_iterator
157 REAL(kind=
dp),
DIMENSION(SIZE(particle_set)) :: at_thread
158 REAL(kind=
dp),
DIMENSION(nexp_max) :: alpha_ppl
159 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: cval_lpot, h1_1block, h1_2block, &
160 h1_3block, h_block, p_block, rpgfa, &
161 rpgfb, sphi_a, sphi_b, zeta, zetb
162 REAL(kind=
dp),
DIMENSION(:),
POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
163 set_radius_a, set_radius_b
164 REAL(kind=
dp),
DIMENSION(4, nexp_max) :: cval_ppl
165 REAL(kind=
dp),
DIMENSION(3, SIZE(particle_set)) :: force_thread
174 do_dr =
PRESENT(deltar)
175 doat =
PRESENT(atcore)
176 IF ((calculate_forces .OR. doat) .AND. do_dr)
THEN
177 cpabort(
"core_ppl: incompatible options")
183 libgrpp_local = .false.
185 IF (calculate_forces)
THEN
186 CALL timeset(routinen//
"_forces", handle)
188 CALL timeset(routinen, handle)
191 nkind =
SIZE(atomic_kind_set)
192 natom =
SIZE(particle_set)
197 IF (
PRESENT(cell_to_index))
THEN
198 cpassert(
ASSOCIATED(cell_to_index))
200 cpabort(
"Missing cell_to_index for k-point calculation")
204 IF (calculate_forces .OR. doat)
THEN
205 IF (
SIZE(matrix_p, 1) == 2)
THEN
207 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
208 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
209 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
210 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
214 force_thread = 0.0_dp
220 maxsgf=maxsgf, maxnset=maxnset, maxlppl=maxlppl, &
221 basis_type=basis_type)
223 maxl = max(maxlgto, maxlppl)
226 ldsab = max(maxco,
ncoset(maxlppl), maxsgf, maxlppl)
227 ldai =
ncoset(maxl + nder + 1)
229 ALLOCATE (basis_set_list(nkind))
231 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
232 IF (
ASSOCIATED(basis_set_a))
THEN
233 basis_set_list(ikind)%gto_basis_set => basis_set_a
235 NULLIFY (basis_set_list(ikind)%gto_basis_set)
239 CALL prepare_upf_local(qs_kind_set, 2*maxlgto + 1, have_upf)
285 ALLOCATE (hab(ldsab, ldsab, maxnset, maxnset), work(ldsab, ldsab*maxder))
287 ALLOCATE (upf_v(ldsab, ldsab))
288 IF (calculate_forces .OR. do_dr)
ALLOCATE (upf_da(ldsab, ldsab, 3), upf_db(ldsab, ldsab, 3))
290 ldai =
ncoset(2*maxlgto + 2*nder)
291 ALLOCATE (ppl_work(ldai, ldai, max(maxder, 2*maxlgto + 2*nder + 1)))
292 IF (calculate_forces .OR. doat)
THEN
293 ALLOCATE (pab(maxco, maxco, maxnset, maxnset))
295 ALLOCATE (ppl_fwork(ldai, ldai, maxder))
299 DO slot = 1, sab_orb(1)%nl_size
302 ALLOCATE (hab2(ldsab, ldsab, 4, maxnset, maxnset))
303 ALLOCATE (hab2_w(ldsab, ldsab, 6))
304 ALLOCATE (ppl_fwork(ldai, ldai, maxder))
307 ikind = sab_orb(1)%nlist_task(slot)%ikind
308 jkind = sab_orb(1)%nlist_task(slot)%jkind
309 iatom = sab_orb(1)%nlist_task(slot)%iatom
310 jatom = sab_orb(1)%nlist_task(slot)%jatom
311 cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
312 rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
314 basis_set_a => basis_set_list(ikind)%gto_basis_set
315 IF (.NOT.
ASSOCIATED(basis_set_a)) cycle
316 basis_set_b => basis_set_list(jkind)%gto_basis_set
317 IF (.NOT.
ASSOCIATED(basis_set_b)) cycle
323 first_sgfa => basis_set_a%first_sgf
324 la_max => basis_set_a%lmax
325 la_min => basis_set_a%lmin
326 npgfa => basis_set_a%npgf
327 nseta = basis_set_a%nset
328 nsgfa => basis_set_a%nsgf_set
329 rpgfa => basis_set_a%pgf_radius
330 set_radius_a => basis_set_a%set_radius
331 sphi_a => basis_set_a%sphi
332 zeta => basis_set_a%zet
334 first_sgfb => basis_set_b%first_sgf
335 lb_max => basis_set_b%lmax
336 lb_min => basis_set_b%lmin
337 npgfb => basis_set_b%npgf
338 nsetb = basis_set_b%nset
339 nsgfb => basis_set_b%nsgf_set
340 rpgfb => basis_set_b%pgf_radius
341 set_radius_b => basis_set_b%set_radius
342 sphi_b => basis_set_b%sphi
343 zetb => basis_set_b%zet
345 dab = sqrt(sum(rab*rab))
348 img = cell_to_index(cellind(1), cellind(2), cellind(3))
354 IF (iatom == jatom)
THEN
361 IF (iatom <= jatom)
THEN
371 NULLIFY (h1_1block, h1_2block, h1_3block)
374 row=irow, col=icol, block=h1_1block, found=found)
376 row=irow, col=icol, block=h1_2block, found=found)
378 row=irow, col=icol, block=h1_3block, found=found)
383 IF (calculate_forces .OR. doat)
THEN
386 IF (
ASSOCIATED(p_block))
THEN
388 ncoa = npgfa(iset)*
ncoset(la_max(iset))
389 sgfa = first_sgfa(1, iset)
391 ncob = npgfb(jset)*
ncoset(lb_max(jset))
392 sgfb = first_sgfb(1, jset)
395 IF (iatom <= jatom)
THEN
396 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
397 p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
399 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
400 transpose(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1)))
403 pab(1:ncoa, 1:ncob, iset, jset) = matmul(work(1:ncoa, 1:nsgfb(jset)), &
404 transpose(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
411 IF (do_dr) hab2 = 0._dp
416 CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
417 sgp_potential=sgp_potential)
418 ecp_semi_local = .false.
420 only_gaussians = .true.
421 IF (
ASSOCIATED(gth_potential))
THEN
423 alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
424 lpot_present=lpotextended, ppl_radius=ppl_radius)
427 nct_ppl(1) =
SIZE(cexp_ppl)
428 cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
429 IF (lpotextended)
THEN
431 nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, &
433 cpassert(nexp_lpot < nexp_max)
434 nexp_ppl = nexp_lpot + 1
435 alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
436 nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
438 cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
441 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
442 CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
443 ppl_radius=ppl_radius)
444 do_upf =
ALLOCATED(sgp_potential%upf_data)
446 cpassert(.NOT. ecp_local .AND. .NOT. ecp_semi_local)
447 ELSE IF (ecp_local)
THEN
448 CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
450 cpassert(nexp_ppl <= nexp_max)
451 nct_ppl(1:nloc) = nrloc(1:nloc)
452 alpha_ppl(1:nloc) = bloc(1:nloc)
453 cval_ppl(1, 1:nloc) = aloc(1:nloc)
454 only_gaussians = .false.
456 CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
458 cpassert(nexp_ppl <= nexp_max)
459 nct_ppl(1:n_local) = 1
460 alpha_ppl(1:n_local) = a_local(1:n_local)
461 cval_ppl(1, 1:n_local) = c_local(1:n_local)
463 IF (ecp_semi_local)
THEN
465 npot=npot, nrpot=nrpot, apot=apot, bpot=bpot)
466 ELSE IF (ecp_local)
THEN
467 IF (sum(abs(aloc(1:nloc))) < 1.0e-12_dp) cycle
479 dac = sqrt(sum(rac*rac))
480 rbc(:) = rac(:) - rab(:)
481 dbc = sqrt(sum(rbc*rbc))
482 IF ((maxval(set_radius_a(:)) + ppl_radius < dac) .OR. &
483 (maxval(set_radius_b(:)) + ppl_radius < dbc))
THEN
488 IF (set_radius_a(iset) + ppl_radius < dac) cycle
489 ncoa = npgfa(iset)*
ncoset(la_max(iset))
490 sgfa = first_sgfa(1, iset)
492 IF (set_radius_b(jset) + ppl_radius < dbc) cycle
493 ncob = npgfb(jset)*
ncoset(lb_max(jset))
494 sgfb = first_sgfb(1, jset)
495 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
498 atk0 = f0*sum(hab(1:ncoa, 1:ncob, iset, jset)* &
499 pab(1:ncoa, 1:ncob, iset, jset))
501 IF (calculate_forces)
THEN
508 la_max(iset), zeta(:npgfa(iset), iset), &
509 lb_max(jset), zetb(:npgfb(jset), jset), rab, rac, &
510 upf_v(:ncoa, :ncob), upf_da(:ncoa, :ncob, :), upf_db(:ncoa, :ncob, :))
511 hab(:ncoa, :ncob, iset, jset) = hab(:ncoa, :ncob, iset, jset) + upf_v(:ncoa, :ncob)
513 force_a(i) = sum(pab(:ncoa, :ncob, iset, jset)*upf_da(:ncoa, :ncob, i))
514 force_b(i) = sum(pab(:ncoa, :ncob, iset, jset)*upf_db(:ncoa, :ncob, i))
516 ELSE IF (only_gaussians)
THEN
518 la_max(iset), la_min(iset), npgfa(iset), &
519 rpgfa(:, iset), zeta(:, iset), &
520 lb_max(jset), lb_min(jset), npgfb(jset), &
521 rpgfb(:, jset), zetb(:, jset), &
522 nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
523 rab, dab, rac, dac, rbc, dbc, &
524 hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
525 force_a, force_b, ppl_fwork)
526 ELSE IF (libgrpp_local)
THEN
529 rpgfa(:, iset), zeta(:, iset), &
530 lb_max(jset), lb_min(jset), npgfb(jset), &
531 rpgfb(:, jset), zetb(:, jset), &
532 nexp_ppl, alpha_ppl, cval_ppl(1, :), nct_ppl, &
533 ppl_radius, rab, dab, rac, dac, dbc, &
534 hab(:, :, iset, jset), pab(:, :, iset, jset), &
539 la_max(iset), la_min(iset), npgfa(iset), &
540 rpgfa(:, iset), zeta(:, iset), &
541 lb_max(jset), lb_min(jset), npgfb(jset), &
542 rpgfb(:, jset), zetb(:, jset), &
543 nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
544 rab, dab, rac, dac, rbc, dbc, &
545 hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
546 force_a, force_b, ppl_fwork)
549 IF (ecp_semi_local)
THEN
553 rpgfa(:, iset), zeta(:, iset), &
554 lb_max(jset), lb_min(jset), npgfb(jset), &
555 rpgfb(:, jset), zetb(:, jset), &
556 slmax, npot, bpot, apot, nrpot, &
557 ppl_radius, rab, dab, rac, dac, dbc, &
558 hab(:, :, iset, jset), pab(:, :, iset, jset), &
566 force_thread(1, iatom) = force_thread(1, iatom) + f0*force_a(1)
567 force_thread(2, iatom) = force_thread(2, iatom) + f0*force_a(2)
568 force_thread(3, iatom) = force_thread(3, iatom) + f0*force_a(3)
569 force_thread(1, katom) = force_thread(1, katom) - f0*force_a(1)
570 force_thread(2, katom) = force_thread(2, katom) - f0*force_a(2)
571 force_thread(3, katom) = force_thread(3, katom) - f0*force_a(3)
573 force_thread(1, jatom) = force_thread(1, jatom) + f0*force_b(1)
574 force_thread(2, jatom) = force_thread(2, jatom) + f0*force_b(2)
575 force_thread(3, jatom) = force_thread(3, jatom) + f0*force_b(3)
576 force_thread(1, katom) = force_thread(1, katom) - f0*force_b(1)
577 force_thread(2, katom) = force_thread(2, katom) - f0*force_b(2)
578 force_thread(3, katom) = force_thread(3, katom) - f0*force_b(3)
587 la_max(iset), zeta(:npgfa(iset), iset), &
588 lb_max(jset), zetb(:npgfb(jset), jset), rab, rac, &
589 upf_v(:ncoa, :ncob), upf_da(:ncoa, :ncob, :), upf_db(:ncoa, :ncob, :))
590 hab(:ncoa, :ncob, iset, jset) = hab(:ncoa, :ncob, iset, jset) + upf_v(:ncoa, :ncob)
592 hab2(:ncoa, :ncob, i, iset, jset) = hab2(:ncoa, :ncob, i, iset, jset) &
593 + upf_da(:ncoa, :ncob, i)*(deltar(i, iatom) - deltar(i, katom)) &
594 + upf_db(:ncoa, :ncob, i)*(deltar(i, jatom) - deltar(i, katom))
599 la_max(iset), la_min(iset), npgfa(iset), &
600 rpgfa(:, iset), zeta(:, iset), &
601 lb_max(jset), lb_min(jset), npgfb(jset), &
602 rpgfb(:, jset), zetb(:, jset), &
603 nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
604 rab, dab, rac, dac, rbc, dbc, &
605 vab=hab(:, :, iset, jset), s=ppl_work, &
606 hab2=hab2(:, :, :, iset, jset), hab2_work=hab2_w, fs=ppl_fwork, &
607 deltar=deltar, iatom=iatom, jatom=jatom, katom=katom)
609 IF (ecp_semi_local)
THEN
611 cpabort(
"Option not implemented")
616 la_max(iset), zeta(:npgfa(iset), iset), &
617 lb_max(jset), zetb(:npgfb(jset), jset), rab, rac, &
619 hab(:ncoa, :ncob, iset, jset) = hab(:ncoa, :ncob, iset, jset) + upf_v(:ncoa, :ncob)
620 ELSE IF (only_gaussians)
THEN
624 la_max(iset), la_min(iset), npgfa(iset), &
625 rpgfa(:, iset), zeta(:, iset), &
626 lb_max(jset), lb_min(jset), npgfb(jset), &
627 rpgfb(:, jset), zetb(:, jset), &
628 nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
629 rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
631 ELSE IF (libgrpp_local)
THEN
635 rpgfa(:, iset), zeta(:, iset), &
636 lb_max(jset), lb_min(jset), npgfb(jset), &
637 rpgfb(:, jset), zetb(:, jset), &
638 nexp_ppl, alpha_ppl, cval_ppl(1, :), nct_ppl, &
639 ppl_radius, rab, dab, rac, dac, dbc, &
640 hab(:, :, iset, jset))
644 la_max(iset), la_min(iset), npgfa(iset), &
645 rpgfa(:, iset), zeta(:, iset), &
646 lb_max(jset), lb_min(jset), npgfb(jset), &
647 rpgfb(:, jset), zetb(:, jset), &
648 nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
649 rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
652 IF (ecp_semi_local)
THEN
656 rpgfa(:, iset), zeta(:, iset), &
657 lb_max(jset), lb_min(jset), npgfb(jset), &
658 rpgfb(:, jset), zetb(:, jset), &
659 slmax, npot, bpot, apot, nrpot, &
660 ppl_radius, rab, dab, rac, dac, dbc, &
661 hab(:, :, iset, jset))
667 atk1 = f0*sum(hab(1:ncoa, 1:ncob, iset, jset)* &
668 pab(1:ncoa, 1:ncob, iset, jset))
669 at_thread(katom) = at_thread(katom) + (atk1 - atk0)
677 IF (.NOT. do_dr)
THEN
679 ncoa = npgfa(iset)*
ncoset(la_max(iset))
680 sgfa = first_sgfa(1, iset)
682 ncob = npgfb(jset)*
ncoset(lb_max(jset))
683 sgfb = first_sgfb(1, jset)
688 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab(1:ncoa, 1:ncob, iset, jset), &
689 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
691 IF (iatom <= jatom)
THEN
692 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
693 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
694 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
696 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
697 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
698 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
706 ncoa = npgfa(iset)*
ncoset(la_max(iset))
707 sgfa = first_sgfa(1, iset)
709 ncob = npgfb(jset)*
ncoset(lb_max(jset))
710 sgfb = first_sgfb(1, jset)
711 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab2(1:ncoa, 1:ncob, 1, iset, jset), &
712 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
715 IF (iatom <= jatom)
THEN
716 h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
717 h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
718 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
721 h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
722 h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
723 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
726 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab2(1:ncoa, 1:ncob, 2, iset, jset), &
727 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
730 IF (iatom <= jatom)
THEN
731 h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
732 h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
733 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
736 h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
737 h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
738 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
741 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab2(1:ncoa, 1:ncob, 3, iset, jset), &
742 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
744 IF (iatom <= jatom)
THEN
745 h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
746 h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
747 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
750 h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
751 h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
752 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
758 IF (do_dr)
DEALLOCATE (hab2, ppl_fwork, hab2_w)
761 DEALLOCATE (hab, work, ppl_work)
762 IF (
ALLOCATED(upf_v))
DEALLOCATE (upf_v)
763 IF (
ALLOCATED(upf_da))
DEALLOCATE (upf_da, upf_db)
764 IF (calculate_forces .OR. doat)
THEN
765 DEALLOCATE (pab, ppl_fwork)
782 DEALLOCATE (basis_set_list)
784 IF (calculate_forces .OR. doat)
THEN
787 IF (
SIZE(matrix_p, 1) == 2)
THEN
789 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
790 alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
791 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
792 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
797 IF (calculate_forces)
THEN
801 atom_a = atom_of_kind(iatom)
802 ikind = kind_of(iatom)
803 force(ikind)%gth_ppl(:, atom_a) = force(ikind)%gth_ppl(:, atom_a) + force_thread(:, iatom)
806 DEALLOCATE (atom_of_kind, kind_of)
809 atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
812 IF (calculate_forces .AND. use_virial)
THEN
813 virial%pv_ppl = virial%pv_ppl + pv_thread
814 virial%pv_virial = virial%pv_virial + pv_thread
817 CALL timestop(handle)
835 qs_kind_set, atomic_kind_set, particle_set, sac_ppl, &
841 LOGICAL,
INTENT(IN) :: calculate_forces
842 LOGICAL :: use_virial
843 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
848 CHARACTER(LEN=*),
INTENT(IN) :: basis_type
850 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_core_ppl_ri'
851 INTEGER,
PARAMETER :: nexp_max = 30
853 INTEGER :: atom_a, handle, i, iatom, ikind, iset, katom, kkind, maxco, maxlgto, maxsgf, &
854 n_local, natom, ncoa, nexp_lpot, nexp_ppl, nfun, nkind, nloc, nseta, sgfa, sgfb, slot
855 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
856 INTEGER,
DIMENSION(1:10) :: nrloc
857 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, nct_lpot, npgfa, nsgfa
858 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa
859 INTEGER,
DIMENSION(nexp_max) :: nct_ppl
860 LOGICAL :: do_upf, ecp_local, ecp_semi_local, &
862 REAL(kind=
dp) :: alpha, dac, ppl_radius
863 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: va, work
864 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dva, dvas
865 REAL(kind=
dp),
DIMENSION(1:10) :: aloc, bloc
866 REAL(kind=
dp),
DIMENSION(3) :: force_a, rac
867 REAL(kind=
dp),
DIMENSION(3, 3) :: pv_thread
871 REAL(kind=
dp),
DIMENSION(nexp_max) :: alpha_ppl
872 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: bcon, cval_lpot, rpgfa, sphi_a, zeta
873 REAL(kind=
dp),
DIMENSION(:),
POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
875 REAL(kind=
dp),
DIMENSION(4, nexp_max) :: cval_ppl
876 REAL(kind=
dp),
DIMENSION(3, SIZE(particle_set)) :: force_thread
884 IF (calculate_forces)
THEN
885 CALL timeset(routinen//
"_forces", handle)
887 CALL timeset(routinen, handle)
890 nkind =
SIZE(atomic_kind_set)
891 natom =
SIZE(particle_set)
893 force_thread = 0.0_dp
897 ALLOCATE (basis_set_list(nkind))
899 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
900 IF (
ASSOCIATED(basis_set))
THEN
901 basis_set_list(ikind)%gto_basis_set => basis_set
903 NULLIFY (basis_set_list(ikind)%gto_basis_set)
907 CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxsgf=maxsgf, maxlgto=maxlgto, basis_type=basis_type)
908 CALL prepare_upf_local(qs_kind_set, maxlgto + 1)
932 ALLOCATE (va(maxco), work(maxsgf))
933 IF (calculate_forces)
THEN
934 ALLOCATE (dva(maxco, 3), dvas(maxco, 3))
938 DO slot = 1, sac_ppl(1)%nl_size
940 ikind = sac_ppl(1)%nlist_task(slot)%ikind
941 kkind = sac_ppl(1)%nlist_task(slot)%jkind
942 iatom = sac_ppl(1)%nlist_task(slot)%iatom
943 katom = sac_ppl(1)%nlist_task(slot)%jatom
944 rac(1:3) = sac_ppl(1)%nlist_task(slot)%r(1:3)
945 atom_a = atom_of_kind(iatom)
947 basis_set => basis_set_list(ikind)%gto_basis_set
948 IF (.NOT.
ASSOCIATED(basis_set)) cycle
951 first_sgfa => basis_set%first_sgf
952 la_max => basis_set%lmax
953 la_min => basis_set%lmin
954 npgfa => basis_set%npgf
955 nseta = basis_set%nset
956 nsgfa => basis_set%nsgf_set
957 nfun = basis_set%nsgf
958 rpgfa => basis_set%pgf_radius
959 set_radius_a => basis_set%set_radius
960 sphi_a => basis_set%sphi
961 zeta => basis_set%zet
963 CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
964 sgp_potential=sgp_potential)
965 ecp_semi_local = .false.
967 IF (
ASSOCIATED(gth_potential))
THEN
969 alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
970 lpot_present=lpotextended, ppl_radius=ppl_radius)
973 nct_ppl(1) =
SIZE(cexp_ppl)
974 cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
975 IF (lpotextended)
THEN
977 nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, cval_lpot=cval_lpot)
978 cpassert(nexp_lpot < nexp_max)
979 nexp_ppl = nexp_lpot + 1
980 alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
981 nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
983 cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
986 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
987 CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
988 ppl_radius=ppl_radius)
989 cpassert(.NOT. ecp_semi_local)
990 do_upf =
ALLOCATED(sgp_potential%upf_data)
992 cpassert(.NOT. ecp_local)
993 ELSE IF (ecp_local)
THEN
994 CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
995 IF (sum(abs(aloc(1:nloc))) < 1.0e-12_dp) cycle
997 cpassert(nexp_ppl <= nexp_max)
998 nct_ppl(1:nloc) = nrloc(1:nloc)
999 alpha_ppl(1:nloc) = bloc(1:nloc)
1000 cval_ppl(1, 1:nloc) = aloc(1:nloc)
1002 CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
1004 cpassert(nexp_ppl <= nexp_max)
1005 nct_ppl(1:n_local) = 1
1006 alpha_ppl(1:n_local) = a_local(1:n_local)
1007 cval_ppl(1, 1:n_local) = c_local(1:n_local)
1013 dac = sqrt(sum(rac*rac))
1014 IF ((maxval(set_radius_a(:)) + ppl_radius < dac)) cycle
1015 IF (calculate_forces) force_a = 0.0_dp
1016 work(1:nfun) = 0.0_dp
1019 IF (set_radius_a(iset) + ppl_radius < dac) cycle
1021 ncoa = npgfa(iset)*
ncoset(la_max(iset))
1023 IF (calculate_forces)
THEN
1024 CALL upf_local_ri(sgp_potential%upf_local, la_max(iset), zeta(:npgfa(iset), iset), &
1025 rac, va(:ncoa), dva(:ncoa, :))
1027 CALL upf_local_ri(sgp_potential%upf_local, la_max(iset), zeta(:npgfa(iset), iset), rac, va(:ncoa))
1029 ELSE IF (calculate_forces)
THEN
1033 la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1034 nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
1039 la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1040 nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
1044 sgfa = first_sgfa(1, iset)
1045 sgfb = sgfa + nsgfa(iset) - 1
1046 ncoa = npgfa(iset)*
ncoset(la_max(iset))
1047 bcon => sphi_a(1:ncoa, sgfa:sgfb)
1048 work(sgfa:sgfb) = matmul(transpose(bcon), va(1:ncoa))
1049 IF (calculate_forces)
THEN
1050 dvas(1:nsgfa(iset), 1:3) = matmul(transpose(bcon), dva(1:ncoa, 1:3))
1051 force_a(1) = force_a(1) + sum(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 1))
1052 force_a(2) = force_a(2) + sum(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 2))
1053 force_a(3) = force_a(3) + sum(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 3))
1058 lri_ppl_coef(ikind)%v_int(atom_a, 1:nfun) = lri_ppl_coef(ikind)%v_int(atom_a, 1:nfun) + work(1:nfun)
1060 IF (calculate_forces)
THEN
1061 force_thread(1, iatom) = force_thread(1, iatom) + force_a(1)
1062 force_thread(2, iatom) = force_thread(2, iatom) + force_a(2)
1063 force_thread(3, iatom) = force_thread(3, iatom) + force_a(3)
1064 force_thread(1, katom) = force_thread(1, katom) - force_a(1)
1065 force_thread(2, katom) = force_thread(2, katom) - force_a(2)
1066 force_thread(3, katom) = force_thread(3, katom) - force_a(3)
1067 IF (use_virial)
THEN
1073 DEALLOCATE (va, work)
1074 IF (calculate_forces)
THEN
1075 DEALLOCATE (dva, dvas)
1080 IF (calculate_forces)
THEN
1082 atom_a = atom_of_kind(iatom)
1083 ikind = kind_of(iatom)
1084 force(ikind)%gth_ppl(1, atom_a) = force(ikind)%gth_ppl(1, atom_a) + force_thread(1, iatom)
1085 force(ikind)%gth_ppl(2, atom_a) = force(ikind)%gth_ppl(2, atom_a) + force_thread(2, iatom)
1086 force(ikind)%gth_ppl(3, atom_a) = force(ikind)%gth_ppl(3, atom_a) + force_thread(3, iatom)
1089 DEALLOCATE (atom_of_kind, kind_of)
1091 IF (calculate_forces .AND. use_virial)
THEN
1092 virial%pv_ppl = virial%pv_ppl + pv_thread
1093 virial%pv_virial = virial%pv_virial + pv_thread
1096 DEALLOCATE (basis_set_list)
1098 CALL timestop(handle)