87 SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
88 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
89 nimages, cell_to_index, basis_type, deltaR, matrix_l, atcore)
91 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_h, matrix_p
94 LOGICAL,
INTENT(IN) :: calculate_forces
101 POINTER :: sab_orb, sap_ppnl
102 REAL(kind=
dp),
INTENT(IN) :: eps_ppnl
103 INTEGER,
INTENT(IN) :: nimages
104 INTEGER,
DIMENSION(:, :, :),
OPTIONAL,
POINTER :: cell_to_index
105 CHARACTER(LEN=*),
INTENT(IN) :: basis_type
106 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
110 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT), &
113 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_core_ppnl'
115 INTEGER :: atom_a, first_col, handle, i, i_dim, iab, iac, iatom, ib, ibc, icol, ikind, &
116 ilist, img, irow, iset, j, jatom, jb, jkind, jneighbor, kac, katom, kbc, kkind, l, &
117 lc_max, lc_min, ldai, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, &
118 maxsgf, na, natom, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, nseta, &
119 nsgfa, prjc, sgfa, slot
120 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
121 INTEGER,
DIMENSION(3) :: cell_b, cell_c
122 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
124 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa
125 LOGICAL :: do_dr, do_gth, do_kp, do_soc, doat, &
127 REAL(kind=
dp) :: atk, dac, f0, ppnl_radius
128 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: radp
129 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: sab, work
130 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ai_work, lab, work_l
131 REAL(kind=
dp),
DIMENSION(1) :: rprjc, zetc
132 REAL(kind=
dp),
DIMENSION(3) :: fa, fb, rab, rac, rbc
133 REAL(kind=
dp),
DIMENSION(3, 3) :: pv_thread
139 TYPE(
alist_type),
POINTER :: alist_ac, alist_bc
140 REAL(kind=
dp),
DIMENSION(SIZE(particle_set)) :: at_thread
141 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: achint, acint, alkint, bchint, bcint, &
143 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: cprj, h_block, l_block_x, l_block_y, &
144 l_block_z, p_block, r_2block, &
145 r_3block, rpgfa, sphi_a, vprj_ppnl, &
147 REAL(kind=
dp),
DIMENSION(:),
POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
148 REAL(kind=
dp),
DIMENSION(3, SIZE(particle_set)) :: force_thread
162 IF (
PRESENT(deltar)) do_dr = .true.
164 IF (
PRESENT(atcore)) doat = .true.
165 IF ((calculate_forces .OR. doat) .AND. do_dr)
THEN
166 cpabort(
"core_ppl: incompatible options")
169 IF (calculate_forces)
THEN
170 CALL timeset(routinen//
"_forces", handle)
172 CALL timeset(routinen, handle)
175 do_soc =
PRESENT(matrix_l)
177 ppnl_present =
ASSOCIATED(sap_ppnl)
179 IF (ppnl_present)
THEN
181 nkind =
SIZE(atomic_kind_set)
182 natom =
SIZE(particle_set)
184 do_kp = (nimages > 1)
187 IF (
PRESENT(cell_to_index))
THEN
188 cpassert(
ASSOCIATED(cell_to_index))
190 cpabort(
"Missing cell_to_index for k-point calculation")
194 IF (calculate_forces .OR. doat)
THEN
195 IF (
SIZE(matrix_p, 1) == 2)
THEN
197 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
198 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
199 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
200 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
213 basis_type=basis_type)
215 maxl = max(maxlgto, maxlppnl)
218 ldsab = max(maxco,
ncoset(maxlppnl), maxsgf, maxppnl)
219 ldai =
ncoset(maxl + nder + 1)
222 ALLOCATE (sap_int(nkind*nkind))
223 DO i = 1, nkind*nkind
224 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
225 sap_int(i)%nalist = 0
229 ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
231 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
232 IF (
ASSOCIATED(orb_basis_set))
THEN
233 basis_set(ikind)%gto_basis_set => orb_basis_set
235 NULLIFY (basis_set(ikind)%gto_basis_set)
237 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
238 NULLIFY (gpotential(ikind)%gth_potential)
239 NULLIFY (spotential(ikind)%sgp_potential)
240 IF (
ASSOCIATED(gth_potential))
THEN
241 gpotential(ikind)%gth_potential => gth_potential
242 IF (do_soc .AND. (.NOT. gth_potential%soc))
THEN
243 cpabort(
"Spin-orbit coupling selected, but GTH potential without SOC parameters provided")
245 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
246 spotential(ikind)%sgp_potential => sgp_potential
251 DO slot = 1, sap_ppnl(1)%nl_size
253 ikind = sap_ppnl(1)%nlist_task(slot)%ikind
254 kkind = sap_ppnl(1)%nlist_task(slot)%jkind
255 iatom = sap_ppnl(1)%nlist_task(slot)%iatom
256 katom = sap_ppnl(1)%nlist_task(slot)%jatom
257 nlist = sap_ppnl(1)%nlist_task(slot)%nlist
258 ilist = sap_ppnl(1)%nlist_task(slot)%ilist
259 nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
261 iac = ikind + nkind*(kkind - 1)
262 IF (.NOT.
ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
263 IF (.NOT.
ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
264 .NOT.
ASSOCIATED(spotential(kkind)%sgp_potential)) cycle
265 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist))
THEN
266 sap_int(iac)%a_kind = ikind
267 sap_int(iac)%p_kind = kkind
268 sap_int(iac)%nalist = nlist
269 ALLOCATE (sap_int(iac)%alist(nlist))
271 NULLIFY (sap_int(iac)%alist(i)%clist)
272 sap_int(iac)%alist(i)%aatom = 0
273 sap_int(iac)%alist(i)%nclist = 0
276 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist(ilist)%clist))
THEN
277 sap_int(iac)%alist(ilist)%aatom = iatom
278 sap_int(iac)%alist(ilist)%nclist = nneighbor
279 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
281 sap_int(iac)%alist(ilist)%clist(i)%catom = 0
299 ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
301 ALLOCATE (ai_work(ldai, ldai,
ncoset(nder + 1)))
304 ALLOCATE (lab(ldsab, ldsab, 3), work_l(ldsab, ldsab, 3))
309 DO slot = 1, sap_ppnl(1)%nl_size
311 ikind = sap_ppnl(1)%nlist_task(slot)%ikind
312 kkind = sap_ppnl(1)%nlist_task(slot)%jkind
313 iatom = sap_ppnl(1)%nlist_task(slot)%iatom
314 katom = sap_ppnl(1)%nlist_task(slot)%jatom
315 nlist = sap_ppnl(1)%nlist_task(slot)%nlist
316 ilist = sap_ppnl(1)%nlist_task(slot)%ilist
317 nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
318 jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
319 cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
320 rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
322 iac = ikind + nkind*(kkind - 1)
323 IF (.NOT.
ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
325 first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
326 la_max => basis_set(ikind)%gto_basis_set%lmax
327 la_min => basis_set(ikind)%gto_basis_set%lmin
328 npgfa => basis_set(ikind)%gto_basis_set%npgf
329 nseta = basis_set(ikind)%gto_basis_set%nset
330 nsgfa = basis_set(ikind)%gto_basis_set%nsgf
331 nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
332 rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
333 set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
334 sphi_a => basis_set(ikind)%gto_basis_set%sphi
335 zeta => basis_set(ikind)%gto_basis_set%zet
337 IF (
ASSOCIATED(gpotential(kkind)%gth_potential))
THEN
340 alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
341 cprj => gpotential(kkind)%gth_potential%cprj
342 lppnl = gpotential(kkind)%gth_potential%lppnl
343 nppnl = gpotential(kkind)%gth_potential%nppnl
344 nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
345 ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
346 vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
347 wprj_ppnl => gpotential(kkind)%gth_potential%wprj_ppnl
348 ELSE IF (
ASSOCIATED(spotential(kkind)%sgp_potential))
THEN
351 nprjc = spotential(kkind)%sgp_potential%nppnl
352 IF (nprjc == 0) cycle
353 nnl = spotential(kkind)%sgp_potential%n_nonlocal
354 lppnl = spotential(kkind)%sgp_potential%lmax
355 a_nl => spotential(kkind)%sgp_potential%a_nonlocal
356 ppnl_radius = spotential(kkind)%sgp_potential%ppnl_radius
358 radp(:) = ppnl_radius
359 cprj => spotential(kkind)%sgp_potential%cprj_ppnl
360 hprj => spotential(kkind)%sgp_potential%vprj_ppnl
361 nppnl =
SIZE(cprj, 2)
366 dac = sqrt(sum(rac*rac))
367 clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
371 ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
372 clist%achint(nsgfa, nppnl, maxder), &
373 clist%alint(nsgfa, nppnl, 3), &
374 clist%alkint(nsgfa, nppnl, 3))
376 clist%achint = 0.0_dp
378 clist%alkint = 0.0_dp
381 NULLIFY (clist%sgf_list)
383 ncoa = npgfa(iset)*
ncoset(la_max(iset))
384 sgfa = first_sgfa(1, iset)
390 nprjc = nprj_ppnl(l)*
nco(l)
391 IF (nprjc == 0) cycle
392 rprjc(1) = ppnl_radius
393 IF (set_radius_a(iset) + rprjc(1) < dac) cycle
394 lc_max = l + 2*(nprj_ppnl(l) - 1)
396 zetc(1) = alpha_ppnl(l)
400 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
401 lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .true., ai_work, ldai)
407 first_col = (i - 1)*ldsab
410 work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
411 matmul(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
417 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
418 lc_max, 1, zetc, rprjc, -rac, [0._dp, 0._dp, 0._dp], lab)
420 work_l(1:na, prjc:prjc + nb - 1, i_dim) = &
421 matmul(lab(1:na, 1:np, i_dim), cprj(1:np, prjc:prjc + nb - 1))
432 first_col = (i - 1)*ldsab + 1
436 clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
437 matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))
441 clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
442 matmul(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))
446 clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
447 matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work_l(1:np, 1:nb, i_dim))
448 clist%alkint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
449 matmul(clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim), wprj_ppnl(1:nb, 1:nb))
455 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
456 lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .true., ai_work, ldai)
461 first_col = (i - 1)*ldsab + 1
465 work(1:np, 1:nb) = matmul(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb))
469 clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
470 matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb))
472 ncoc = sgfa + nsgf_seta(iset) - 1
474 clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
479 clist%maxac = maxval(abs(clist%acint(:, :, 1)))
480 clist%maxach = maxval(abs(clist%achint(:, :, 1)))
481 IF (.NOT. do_gth)
DEALLOCATE (radp)
484 DEALLOCATE (sab, ai_work, work)
485 IF (do_soc)
DEALLOCATE (lab, work_l)
493 force_thread = 0.0_dp
524 DO slot = 1, sab_orb(1)%nl_size
526 ikind = sab_orb(1)%nlist_task(slot)%ikind
527 jkind = sab_orb(1)%nlist_task(slot)%jkind
528 iatom = sab_orb(1)%nlist_task(slot)%iatom
529 jatom = sab_orb(1)%nlist_task(slot)%jatom
530 cell_b(:) = sab_orb(1)%nlist_task(slot)%cell(:)
531 rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
533 IF (.NOT.
ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
534 IF (.NOT.
ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
536 iab = ikind + nkind*(jkind - 1)
539 IF (iatom == jatom)
THEN
546 img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
552 IF (iatom <= jatom)
THEN
562 NULLIFY (l_block_x, l_block_y, l_block_z)
569 NULLIFY (r_2block, r_3block)
574 IF (calculate_forces .OR. doat)
THEN
580 IF (
ASSOCIATED(h_block))
THEN
585 iac = ikind + nkind*(kkind - 1)
586 ibc = jkind + nkind*(kkind - 1)
587 IF (.NOT.
ASSOCIATED(sap_int(iac)%alist)) cycle
588 IF (.NOT.
ASSOCIATED(sap_int(ibc)%alist)) cycle
589 CALL get_alist(sap_int(iac), alist_ac, iatom)
590 CALL get_alist(sap_int(ibc), alist_bc, jatom)
591 IF (.NOT.
ASSOCIATED(alist_ac)) cycle
592 IF (.NOT.
ASSOCIATED(alist_bc)) cycle
593 DO kac = 1, alist_ac%nclist
594 DO kbc = 1, alist_bc%nclist
595 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
596 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0))
THEN
597 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
598 acint => alist_ac%clist(kac)%acint
599 bcint => alist_bc%clist(kbc)%acint
600 achint => alist_ac%clist(kac)%achint
601 bchint => alist_bc%clist(kbc)%achint
603 alkint => alist_ac%clist(kac)%alkint
604 blkint => alist_bc%clist(kbc)%alkint
610 IF (.NOT. do_dr)
THEN
611 IF (iatom <= jatom)
THEN
612 h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
613 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 1)))
615 h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
616 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 1)))
620 IF (iatom <= jatom)
THEN
621 l_block_x(1:na, 1:nb) = l_block_x(1:na, 1:nb) + &
622 matmul(alkint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 1)))
623 l_block_y(1:na, 1:nb) = l_block_y(1:na, 1:nb) + &
624 matmul(alkint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))
625 l_block_z(1:na, 1:nb) = l_block_z(1:na, 1:nb) + &
626 matmul(alkint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))
629 l_block_x(1:nb, 1:na) = l_block_x(1:nb, 1:na) + &
630 matmul(blkint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 1)))
631 l_block_y(1:nb, 1:na) = l_block_y(1:nb, 1:na) + &
632 matmul(blkint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))
633 l_block_z(1:nb, 1:na) = l_block_z(1:nb, 1:na) + &
634 matmul(blkint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))
638 IF (calculate_forces)
THEN
639 IF (
ASSOCIATED(p_block))
THEN
640 katom = alist_ac%clist(kac)%catom
643 IF (iatom <= jatom)
THEN
644 fa(i) = sum(p_block(1:na, 1:nb)* &
645 matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1))))
646 fb(i) = sum(p_block(1:na, 1:nb)* &
647 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j))))
649 fa(i) = sum(p_block(1:nb, 1:na)* &
650 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j))))
651 fb(i) = sum(p_block(1:nb, 1:na)* &
652 matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1))))
654 force_thread(i, iatom) = force_thread(i, iatom) + f0*fa(i)
655 force_thread(i, katom) = force_thread(i, katom) - f0*fa(i)
656 force_thread(i, jatom) = force_thread(i, jatom) + f0*fb(i)
657 force_thread(i, katom) = force_thread(i, katom) - f0*fb(i)
661 rac = alist_ac%clist(kac)%rac
662 rbc = alist_bc%clist(kbc)%rac
671 katom = alist_ac%clist(kac)%catom
672 IF (iatom <= jatom)
THEN
673 h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
674 (deltar(i, iatom) - deltar(i, katom))* &
675 matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1)))
677 h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
678 (deltar(i, jatom) - deltar(i, katom))* &
679 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j)))
681 h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
682 (deltar(i, iatom) - deltar(i, katom))* &
683 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j)))
684 h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
685 (deltar(i, jatom) - deltar(i, katom))* &
686 matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1)))
690 katom = alist_ac%clist(kac)%catom
691 IF (iatom <= jatom)
THEN
692 r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
693 (deltar(i, iatom) - deltar(i, katom))* &
694 matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1)))
696 r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
697 (deltar(i, jatom) - deltar(i, katom))* &
698 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j)))
700 r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
701 (deltar(i, iatom) - deltar(i, katom))* &
702 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j)))
703 r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
704 (deltar(i, jatom) - deltar(i, katom))* &
705 matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1)))
709 katom = alist_ac%clist(kac)%catom
710 IF (iatom <= jatom)
THEN
711 r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
712 (deltar(i, iatom) - deltar(i, katom))* &
713 matmul(acint(1:na, 1:np, j), transpose(bchint(1:nb, 1:np, 1)))
715 r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
716 (deltar(i, jatom) - deltar(i, katom))* &
717 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, j)))
719 r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
720 (deltar(i, iatom) - deltar(i, katom))* &
721 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, j)))
722 r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
723 (deltar(i, jatom) - deltar(i, katom))* &
724 matmul(bcint(1:nb, 1:np, j), transpose(achint(1:na, 1:np, 1)))
729 IF (
ASSOCIATED(p_block))
THEN
730 katom = alist_ac%clist(kac)%catom
731 IF (iatom <= jatom)
THEN
732 atk = sum(p_block(1:na, 1:nb)* &
733 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 1))))
735 atk = sum(p_block(1:nb, 1:na)* &
736 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 1))))
738 at_thread(katom) = at_thread(katom) + f0*atk
763 DEALLOCATE (basis_set, gpotential, spotential)
764 IF (calculate_forces)
THEN
768 atom_a = atom_of_kind(iatom)
769 ikind = kind_of(iatom)
770 force(ikind)%gth_ppnl(:, atom_a) = force(ikind)%gth_ppnl(:, atom_a) + force_thread(:, iatom)
773 DEALLOCATE (atom_of_kind, kind_of)
776 IF (calculate_forces .AND. use_virial)
THEN
777 virial%pv_ppnl = virial%pv_ppnl + pv_thread
778 virial%pv_virial = virial%pv_virial + pv_thread
782 atcore(1:natom) = atcore(1:natom) + at_thread
785 IF (calculate_forces .OR. doat)
THEN
788 IF (
SIZE(matrix_p, 1) == 2)
THEN
790 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
791 alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
792 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
793 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
800 CALL timestop(handle)