(git:6344619)
Loading...
Searching...
No Matches
sap_kind_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief General overlap type integrals containers
9!> \par History
10!> - rewrite of PPNL and OCE integrals
11! **************************************************************************************************
12MODULE sap_kind_types
13
14 USE ai_moments, ONLY: moment
15 USE ai_overlap, ONLY: overlap
16 USE basis_set_types, ONLY: gto_basis_set_p_type, &
17 gto_basis_set_type
18 USE cell_types, ONLY: cell_type, &
19 pbc
20 USE external_potential_types, ONLY: gth_potential_p_type, &
21 gth_potential_type, &
22 sgp_potential_p_type, &
23 sgp_potential_type
24 USE kinds, ONLY: dp
25 USE orbital_pointers, ONLY: init_orbital_pointers, &
26 nco, &
27 ncoset
28 USE particle_types, ONLY: particle_type
29 USE qs_kind_types, ONLY: get_qs_kind, &
30 get_qs_kind_set, &
31 qs_kind_type
32 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
33 USE upf_projector_integrals, ONLY: upf_projector_init, &
34 upf_projector_overlap, &
35 upf_projector_type
36 USE util, ONLY: locate, &
37 sort
38#include "./base/base_uses.f90"
39
40 IMPLICIT NONE
41
42 PRIVATE
43
44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'sap_kind_types'
45
46 TYPE clist_type
47 INTEGER :: catom = -1
48 INTEGER :: nsgf_cnt = -1
49 INTEGER, DIMENSION(:), POINTER :: sgf_list => null()
50 INTEGER, DIMENSION(3) :: cell = -1
51 LOGICAL :: sgf_soft_only = .false.
52 REAL(kind=dp) :: maxac = 0.0_dp
53 REAL(kind=dp) :: maxach = 0.0_dp
54 REAL(kind=dp), DIMENSION(3) :: rac = 0.0_dp
55 REAL (kind=dp), DIMENSION(:, :, :), POINTER :: acint => null()
56 REAL (kind=dp), DIMENSION(:, :, :), POINTER :: achint => null()
57 REAL (kind=dp), DIMENSION(:, :, :), POINTER :: alint => null()
58 REAL (kind=dp), DIMENSION(:, :, :), POINTER :: alkint => null()
59 END TYPE clist_type
60
61 TYPE alist_type
62 INTEGER :: aatom = -1
63 INTEGER :: nclist = -1
64 TYPE(clist_type), DIMENSION(:), POINTER :: clist => null()
65 END TYPE alist_type
66
67 TYPE sap_int_type
68 INTEGER :: a_kind = -1
69 INTEGER :: p_kind = -1
70 INTEGER :: nalist = -1
71 TYPE(alist_type), DIMENSION(:), POINTER :: alist => null()
72 INTEGER, DIMENSION(:), POINTER :: asort => null()
73 INTEGER, DIMENSION(:), POINTER :: aindex => null()
74 END TYPE sap_int_type
75 TYPE clist_complex_type
76 INTEGER :: catom = -1
77 INTEGER :: nsgf_cnt = -1
78 INTEGER, DIMENSION(:), POINTER :: sgf_list => null()
79 INTEGER, DIMENSION(3) :: cell = -1
80 LOGICAL :: sgf_soft_only = .false.
81 REAL(kind=dp) :: maxac = 0.0_dp
82 REAL(kind=dp) :: maxach = 0.0_dp
83 REAL(kind=dp), DIMENSION(3) :: rac = 0.0_dp
84 COMPLEX (KIND=dp), DIMENSION(:, :, :), POINTER :: acint => null()
85 COMPLEX (KIND=dp), DIMENSION(:, :, :), POINTER :: achint => null()
86 COMPLEX (KIND=dp), DIMENSION(:, :, :), POINTER :: alint => null()
87 COMPLEX (KIND=dp), DIMENSION(:, :, :), POINTER :: alkint => null()
88 END TYPE clist_complex_type
89
90 TYPE alist_complex_type
91 INTEGER :: aatom = -1
92 INTEGER :: nclist = -1
93 TYPE(clist_complex_type), DIMENSION(:), POINTER :: clist => null()
94 END TYPE alist_complex_type
95
96 TYPE sap_int_complex_type
97 INTEGER :: a_kind = -1
98 INTEGER :: p_kind = -1
99 INTEGER :: nalist = -1
100 TYPE(alist_complex_type), DIMENSION(:), POINTER :: alist => null()
101 INTEGER, DIMENSION(:), POINTER :: asort => null()
102 INTEGER, DIMENSION(:), POINTER :: aindex => null()
103 END TYPE sap_int_complex_type
104
105 PUBLIC :: sap_int_type, clist_type, alist_type, &
106 sap_int_complex_type, clist_complex_type, alist_complex_type, &
107 release_sap_int, get_alist, alist_pre_align_blk, &
108 alist_post_align_blk, sap_sort, build_sap_ints
109
110 INTERFACE release_sap_int
111 MODULE PROCEDURE release_sap_int_real, release_sap_int_complex
112 END INTERFACE release_sap_int
113 INTERFACE get_alist
114 MODULE PROCEDURE get_alist_real, get_alist_complex
115 END INTERFACE get_alist
116 INTERFACE sap_sort
117 MODULE PROCEDURE sap_sort_real, sap_sort_complex
118 END INTERFACE sap_sort
119
120CONTAINS
121
122!==========================================================================================================
123
124! **************************************************************************************************
125!> \brief ...
126!> \param sap_int ...
127! **************************************************************************************************
128 SUBROUTINE release_sap_int_real (sap_int)
129
130 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
131
132 INTEGER :: i, j, k
133 TYPE(clist_type), POINTER :: clist
134
135 cpassert(ASSOCIATED(sap_int))
136
137 DO i = 1, SIZE(sap_int)
138 IF (ASSOCIATED(sap_int(i)%alist)) THEN
139 DO j = 1, SIZE(sap_int(i)%alist)
140 IF (ASSOCIATED(sap_int(i)%alist(j)%clist)) THEN
141 DO k = 1, SIZE(sap_int(i)%alist(j)%clist)
142 clist => sap_int(i)%alist(j)%clist(k)
143 IF (ASSOCIATED(clist%acint)) THEN
144 DEALLOCATE (clist%acint)
145 END IF
146 IF (ASSOCIATED(clist%sgf_list)) THEN
147 DEALLOCATE (clist%sgf_list)
148 END IF
149 IF (ASSOCIATED(clist%achint)) THEN
150 DEALLOCATE (clist%achint)
151 END IF
152 IF (ASSOCIATED(clist%alint)) THEN
153 DEALLOCATE (clist%alint)
154 END IF
155 IF (ASSOCIATED(clist%alkint)) THEN
156 DEALLOCATE (clist%alkint)
157 END IF
158 END DO
159 DEALLOCATE (sap_int(i)%alist(j)%clist)
160 END IF
161 END DO
162 DEALLOCATE (sap_int(i)%alist)
163 END IF
164 IF (ASSOCIATED(sap_int(i)%asort)) THEN
165 DEALLOCATE (sap_int(i)%asort)
166 END IF
167 IF (ASSOCIATED(sap_int(i)%aindex)) THEN
168 DEALLOCATE (sap_int(i)%aindex)
169 END IF
170 END DO
171
172 DEALLOCATE (sap_int)
173
174 END SUBROUTINE release_sap_int_real
175! **************************************************************************************************
176!> \brief ...
177!> \param sap_int ...
178! **************************************************************************************************
179 SUBROUTINE release_sap_int_complex (sap_int)
180
181 TYPE(sap_int_complex_type), DIMENSION(:), POINTER :: sap_int
182
183 INTEGER :: i, j, k
184 TYPE(clist_complex_type), POINTER :: clist
185
186 cpassert(ASSOCIATED(sap_int))
187
188 DO i = 1, SIZE(sap_int)
189 IF (ASSOCIATED(sap_int(i)%alist)) THEN
190 DO j = 1, SIZE(sap_int(i)%alist)
191 IF (ASSOCIATED(sap_int(i)%alist(j)%clist)) THEN
192 DO k = 1, SIZE(sap_int(i)%alist(j)%clist)
193 clist => sap_int(i)%alist(j)%clist(k)
194 IF (ASSOCIATED(clist%acint)) THEN
195 DEALLOCATE (clist%acint)
196 END IF
197 IF (ASSOCIATED(clist%sgf_list)) THEN
198 DEALLOCATE (clist%sgf_list)
199 END IF
200 IF (ASSOCIATED(clist%achint)) THEN
201 DEALLOCATE (clist%achint)
202 END IF
203 IF (ASSOCIATED(clist%alint)) THEN
204 DEALLOCATE (clist%alint)
205 END IF
206 IF (ASSOCIATED(clist%alkint)) THEN
207 DEALLOCATE (clist%alkint)
208 END IF
209 END DO
210 DEALLOCATE (sap_int(i)%alist(j)%clist)
211 END IF
212 END DO
213 DEALLOCATE (sap_int(i)%alist)
214 END IF
215 IF (ASSOCIATED(sap_int(i)%asort)) THEN
216 DEALLOCATE (sap_int(i)%asort)
217 END IF
218 IF (ASSOCIATED(sap_int(i)%aindex)) THEN
219 DEALLOCATE (sap_int(i)%aindex)
220 END IF
221 END DO
222
223 DEALLOCATE (sap_int)
224
225 END SUBROUTINE release_sap_int_complex
226
227! **************************************************************************************************
228!> \brief ...
229!> \param sap_int ...
230!> \param alist ...
231!> \param atom ...
232! **************************************************************************************************
233 SUBROUTINE get_alist_real (sap_int, alist, atom)
234
235 TYPE(sap_int_type), INTENT(IN) :: sap_int
236 TYPE(alist_type), INTENT(OUT), POINTER :: alist
237 INTEGER, INTENT(IN) :: atom
238
239 INTEGER :: i
240
241 NULLIFY (alist)
242 i = locate(sap_int%asort, atom)
243 IF (i > 0 .AND. i <= SIZE(sap_int%alist)) THEN
244 i = sap_int%aindex(i)
245 alist => sap_int%alist(i)
246 ELSE IF (i == 0) THEN
247 NULLIFY (alist)
248 ELSE
249 cpabort("Invalid result of locate(sap_int%asort, atom)")
250 END IF
251
252 END SUBROUTINE get_alist_real
253! **************************************************************************************************
254!> \brief ...
255!> \param sap_int ...
256!> \param alist ...
257!> \param atom ...
258! **************************************************************************************************
259 SUBROUTINE get_alist_complex (sap_int, alist, atom)
260
261 TYPE(sap_int_complex_type), INTENT(IN) :: sap_int
262 TYPE(alist_complex_type), INTENT(OUT), POINTER :: alist
263 INTEGER, INTENT(IN) :: atom
264
265 INTEGER :: i
266
267 NULLIFY (alist)
268 i = locate(sap_int%asort, atom)
269 IF (i > 0 .AND. i <= SIZE(sap_int%alist)) THEN
270 i = sap_int%aindex(i)
271 alist => sap_int%alist(i)
272 ELSE IF (i == 0) THEN
273 NULLIFY (alist)
274 ELSE
275 cpabort("Invalid result of locate(sap_int%asort, atom)")
276 END IF
277
278 END SUBROUTINE get_alist_complex
279
280! **************************************************************************************************
281!> \brief ...
282!> \param blk_in ...
283!> \param ldin ...
284!> \param blk_out ...
285!> \param ldout ...
286!> \param ilist ...
287!> \param in ...
288!> \param jlist ...
289!> \param jn ...
290! **************************************************************************************************
291 SUBROUTINE alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
292 INTEGER, INTENT(IN) :: in, ilist(*), ldout
293 REAL(dp), INTENT(INOUT) :: blk_out(ldout, *)
294 INTEGER, INTENT(IN) :: ldin
295 REAL(dp), INTENT(IN) :: blk_in(ldin, *)
296 INTEGER, INTENT(IN) :: jlist(*), jn
297
298 INTEGER :: i, i0, i1, i2, i3, inn, inn1, j, j0
299
300 inn = mod(in, 4)
301 inn1 = inn + 1
302 DO j = 1, jn
303 j0 = jlist(j)
304 DO i = 1, inn
305 i0 = ilist(i)
306 blk_out(i, j) = blk_in(i0, j0)
307 END DO
308 DO i = inn1, in, 4
309 i0 = ilist(i)
310 i1 = ilist(i + 1)
311 i2 = ilist(i + 2)
312 i3 = ilist(i + 3)
313 blk_out(i, j) = blk_in(i0, j0)
314 blk_out(i + 1, j) = blk_in(i1, j0)
315 blk_out(i + 2, j) = blk_in(i2, j0)
316 blk_out(i + 3, j) = blk_in(i3, j0)
317 END DO
318 END DO
319 END SUBROUTINE alist_pre_align_blk
320
321! **************************************************************************************************
322!> \brief ...
323!> \param blk_in ...
324!> \param ldin ...
325!> \param blk_out ...
326!> \param ldout ...
327!> \param ilist ...
328!> \param in ...
329!> \param jlist ...
330!> \param jn ...
331! **************************************************************************************************
332 SUBROUTINE alist_post_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
333 INTEGER, INTENT(IN) :: in, ilist(*), ldout
334 REAL(dp), INTENT(INOUT) :: blk_out(ldout, *)
335 INTEGER, INTENT(IN) :: ldin
336 REAL(dp), INTENT(IN) :: blk_in(ldin, *)
337 INTEGER, INTENT(IN) :: jlist(*), jn
338
339 INTEGER :: i, i0, i1, i2, i3, inn, inn1, j, j0
340
341 inn = mod(in, 4)
342 inn1 = inn + 1
343 DO j = 1, jn
344 j0 = jlist(j)
345 DO i = 1, inn
346 i0 = ilist(i)
347 blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j)
348 END DO
349 DO i = inn1, in, 4
350 i0 = ilist(i)
351 i1 = ilist(i + 1)
352 i2 = ilist(i + 2)
353 i3 = ilist(i + 3)
354 blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j)
355 blk_out(i1, j0) = blk_out(i1, j0) + blk_in(i + 1, j)
356 blk_out(i2, j0) = blk_out(i2, j0) + blk_in(i + 2, j)
357 blk_out(i3, j0) = blk_out(i3, j0) + blk_in(i + 3, j)
358 END DO
359 END DO
360 END SUBROUTINE alist_post_align_blk
361
362! **************************************************************************************************
363!> \brief ...
364!> \param sap_int ...
365! **************************************************************************************************
366 SUBROUTINE sap_sort_real (sap_int)
367 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
368
369 INTEGER :: iac, na
370
371! *** Set up a sorting index
372
373!$OMP PARALLEL DEFAULT(NONE) SHARED(sap_int) PRIVATE(iac,na)
374!$OMP DO
375 DO iac = 1, SIZE(sap_int)
376 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
377 na = SIZE(sap_int(iac)%alist)
378 ALLOCATE (sap_int(iac)%asort(na), sap_int(iac)%aindex(na))
379 sap_int(iac)%asort(1:na) = sap_int(iac)%alist(1:na)%aatom
380 CALL sort(sap_int(iac)%asort, na, sap_int(iac)%aindex)
381 END DO
382!$OMP END PARALLEL
383
384 END SUBROUTINE sap_sort_real
385! **************************************************************************************************
386!> \brief ...
387!> \param sap_int ...
388! **************************************************************************************************
389 SUBROUTINE sap_sort_complex (sap_int)
390 TYPE(sap_int_complex_type), DIMENSION(:), POINTER :: sap_int
391
392 INTEGER :: iac, na
393
394! *** Set up a sorting index
395
396!$OMP PARALLEL DEFAULT(NONE) SHARED(sap_int) PRIVATE(iac,na)
397!$OMP DO
398 DO iac = 1, SIZE(sap_int)
399 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
400 na = SIZE(sap_int(iac)%alist)
401 ALLOCATE (sap_int(iac)%asort(na), sap_int(iac)%aindex(na))
402 sap_int(iac)%asort(1:na) = sap_int(iac)%alist(1:na)%aatom
403 CALL sort(sap_int(iac)%asort, na, sap_int(iac)%aindex)
404 END DO
405!$OMP END PARALLEL
406
407 END SUBROUTINE sap_sort_complex
408
409!==========================================================================================================
410
411! **************************************************************************************************
412!> \brief Calculate overlap and optionally momenta <a|x^n|p> between GTOs and nl. pseudo potential projectors
413!> adapted from core_ppnl.F::build_core_ppnl
414!> \param sap_int allocated in parent routine (nkind*nkind)
415!> \param sap_ppnl ...
416!> \param qs_kind_set ...
417!> \param nder Either number of derivatives or order of moments
418!> \param moment_mode if present and true, moments are calculated instead of derivatives
419!> \param refpoint optionally the reference point for moment calculation
420!> \param particle_set needed if refpoint is present
421!> \param cell needed if refpoint is present
422!> \param pseudoatom If we want to constrain the calculations to katom == pseudoatom
423! **************************************************************************************************
424 SUBROUTINE build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder, moment_mode, refpoint, particle_set, cell, pseudoatom)
425 TYPE(sap_int_type), DIMENSION(:), INTENT(INOUT), &
426 POINTER :: sap_int
427 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
428 INTENT(IN), POINTER :: sap_ppnl
429 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
430 POINTER :: qs_kind_set
431 INTEGER, INTENT(IN) :: nder
432 LOGICAL, INTENT(IN), OPTIONAL :: moment_mode
433 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
434 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
435 OPTIONAL, POINTER :: particle_set
436 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
437 INTEGER, OPTIONAL :: pseudoatom
438
439 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_sap_ints'
440
441 INTEGER :: first_col, handle, i, iac, iatom, ikind, ilist, iset, j, jneighbor, katom, kkind, &
442 l, lc_max, lc_min, ldai, ldai_nl, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, &
443 maxppnl, maxsgf, na, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, &
444 nseta, nsgfa, prjc, sgfa, slot
445 INTEGER, DIMENSION(3) :: cell_c
446 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
447 nsgf_seta
448 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
449 LOGICAL :: do_upf, dogth, my_momentmode, my_ref
450 LOGICAL, DIMENSION(0:9) :: is_nonlocal
451 REAL(kind=dp) :: dac, ppnl_radius
452 REAL(kind=dp), ALLOCATABLE :: upf_int(:, :, :)
453 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: radp
454 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
455 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work
456 REAL(kind=dp), DIMENSION(1) :: rprjc, zetc
457 REAL(kind=dp), DIMENSION(3) :: rac, raf, rcf, rf
458 REAL(kind=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
459 REAL(kind=dp), DIMENSION(:, :), POINTER :: cprj, h_nl, rpgfa, sphi_a, vprj_ppnl, &
460 zeta
461 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: c_nl
462 TYPE(clist_type), POINTER :: clist
463 TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
464 TYPE(gth_potential_type), POINTER :: gth_potential
465 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
466 DIMENSION(:) :: basis_set
467 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
468 TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
469 TYPE(sgp_potential_type), POINTER :: sgp_potential
470 TYPE(upf_projector_type), POINTER :: upf
471
472 CALL timeset(routinen, handle)
473
474 nkind = SIZE(qs_kind_set)
475 maxder = ncoset(nder)
476
477 ! determine whether moments or derivatives should be calculated
478 my_momentmode = .false.
479 IF (PRESENT(moment_mode)) THEN
480 my_momentmode = moment_mode
481 END IF
482
483 my_ref = .false.
484 IF (PRESENT(refpoint)) THEN
485 ! A reference point also requires the cell and particle coordinates.
486 cpassert((PRESENT(cell) .AND. PRESENT(particle_set)))
487 rf = refpoint
488 my_ref = .true.
489 END IF
490
491 CALL get_qs_kind_set(qs_kind_set, &
492 maxco=maxco, &
493 maxlgto=maxlgto, &
494 maxsgf=maxsgf, &
495 maxlppnl=maxlppnl, &
496 maxppnl=maxppnl)
497
498 maxl = max(maxlgto, maxlppnl)
499 CALL init_orbital_pointers(maxl + nder + 1)
500
501 ldsab = max(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
502 IF (.NOT. my_momentmode) THEN
503 ldai = ncoset(maxl + nder + 1)
504 ELSE
505 ldai = maxco
506 END IF
507
508 !set up direct access to basis and potential
509 ldai_nl = 0
510 NULLIFY (gpotential, spotential)
511 ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
512 DO ikind = 1, nkind
513 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
514 IF (ASSOCIATED(orb_basis_set)) THEN
515 basis_set(ikind)%gto_basis_set => orb_basis_set
516 ELSE
517 NULLIFY (basis_set(ikind)%gto_basis_set)
518 END IF
519 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
520 NULLIFY (gpotential(ikind)%gth_potential)
521 NULLIFY (spotential(ikind)%sgp_potential)
522 IF (ASSOCIATED(gth_potential)) THEN
523 gpotential(ikind)%gth_potential => gth_potential
524 ldai_nl = max(ldai_nl, ncoset(gth_potential%lprj_ppnl_max))
525 ELSE IF (ASSOCIATED(sgp_potential)) THEN
526 spotential(ikind)%sgp_potential => sgp_potential
527 IF (ALLOCATED(sgp_potential%upf_data)) THEN
528 IF (sgp_potential%upf_data%number_of_proj > 0) THEN
529 IF (sgp_potential%upf_projectors%degree < maxlgto + nder) THEN
530 CALL upf_projector_init(sgp_potential%upf_data, maxlgto + nder, sgp_potential%upf_projectors)
531 END IF
532 END IF
533 END IF
534 ldai_nl = max(ldai_nl, sgp_potential%n_nonlocal*ncoset(sgp_potential%lmax))
535 END IF
536 END DO
537
538 !allocate sap int
539 DO slot = 1, sap_ppnl(1)%nl_size
540
541 ikind = sap_ppnl(1)%nlist_task(slot)%ikind
542 kkind = sap_ppnl(1)%nlist_task(slot)%jkind
543 iatom = sap_ppnl(1)%nlist_task(slot)%iatom
544 katom = sap_ppnl(1)%nlist_task(slot)%jatom
545 nlist = sap_ppnl(1)%nlist_task(slot)%nlist
546 ilist = sap_ppnl(1)%nlist_task(slot)%ilist
547 nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
548
549 iac = ikind + nkind*(kkind - 1)
550 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
551 IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
552 .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) cycle
553 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
554 sap_int(iac)%a_kind = ikind
555 sap_int(iac)%p_kind = kkind
556 sap_int(iac)%nalist = nlist
557 ALLOCATE (sap_int(iac)%alist(nlist))
558 DO i = 1, nlist
559 NULLIFY (sap_int(iac)%alist(i)%clist)
560 sap_int(iac)%alist(i)%aatom = 0
561 sap_int(iac)%alist(i)%nclist = 0
562 END DO
563 END IF
564 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
565 sap_int(iac)%alist(ilist)%aatom = iatom
566 sap_int(iac)%alist(ilist)%nclist = nneighbor
567 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
568 DO i = 1, nneighbor
569 sap_int(iac)%alist(ilist)%clist(i)%catom = 0
570 END DO
571 END IF
572 END DO
573
574 !calculate the overlap integrals <a|pp>
575!$OMP PARALLEL &
576!$OMP DEFAULT (NONE) &
577!$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, my_momentmode, ldai_nl, &
578!$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco, my_ref, cell, particle_set, rf, pseudoatom) &
579!$OMP PRIVATE (ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
580!$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
581!$OMP slot, sphi_a, zeta, cprj, hprj, lppnl, nppnl, nprj_ppnl, &
582!$OMP do_upf, upf, upf_int, &
583!$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, &
584!$OMP ncoc, rpgfa, first_col, vprj_ppnl, i, j, l, dogth, &
585!$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl, &
586!$OMP na, nb, np, nnl, is_nonlocal, a_nl, h_nl, c_nl, radp, raf, rcf)
587
588 ! allocate work storage
589 ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
590 sab = 0.0_dp
591 IF (.NOT. my_momentmode) THEN
592 ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
593 ELSE
594 ALLOCATE (ai_work(ldai, ldai_nl, ncoset(nder + 1)))
595 END IF
596 ai_work = 0.0_dp
597
598!$OMP DO SCHEDULE(GUIDED)
599 ! loop over neighbourlist
600 DO slot = 1, sap_ppnl(1)%nl_size
601 ikind = sap_ppnl(1)%nlist_task(slot)%ikind
602 kkind = sap_ppnl(1)%nlist_task(slot)%jkind
603 iatom = sap_ppnl(1)%nlist_task(slot)%iatom
604 katom = sap_ppnl(1)%nlist_task(slot)%jatom
605 nlist = sap_ppnl(1)%nlist_task(slot)%nlist
606 ilist = sap_ppnl(1)%nlist_task(slot)%ilist
607 nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
608 jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
609 cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
610 rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
611
612 iac = ikind + nkind*(kkind - 1)
613 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
614 ! get definition of basis set
615 first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
616 la_max => basis_set(ikind)%gto_basis_set%lmax
617 la_min => basis_set(ikind)%gto_basis_set%lmin
618 npgfa => basis_set(ikind)%gto_basis_set%npgf
619 nseta = basis_set(ikind)%gto_basis_set%nset
620 nsgfa = basis_set(ikind)%gto_basis_set%nsgf
621 nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
622 rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
623 set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
624 sphi_a => basis_set(ikind)%gto_basis_set%sphi
625 zeta => basis_set(ikind)%gto_basis_set%zet
626 ! get definition of PP projectors
627 do_upf = .false.
628 IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN
629 ! GTH potential
630 dogth = .true.
631 alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
632 cprj => gpotential(kkind)%gth_potential%cprj
633 lppnl = gpotential(kkind)%gth_potential%lppnl
634 nppnl = gpotential(kkind)%gth_potential%nppnl
635 nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
636 ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
637 vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
638 ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN
639 ! SGP potential
640 dogth = .false.
641 nprjc = spotential(kkind)%sgp_potential%nppnl
642 IF (nprjc == 0) cycle
643 do_upf = ALLOCATED(spotential(kkind)%sgp_potential%upf_data)
644 IF (do_upf) THEN
645 upf => spotential(kkind)%sgp_potential%upf_projectors
646 nppnl = upf%nproj
647 ELSE
648 nnl = spotential(kkind)%sgp_potential%n_nonlocal
649 lppnl = spotential(kkind)%sgp_potential%lmax
650 is_nonlocal = .false.
651 is_nonlocal(0:lppnl) = spotential(kkind)%sgp_potential%is_nonlocal(0:lppnl)
652 a_nl => spotential(kkind)%sgp_potential%a_nonlocal
653 h_nl => spotential(kkind)%sgp_potential%h_nonlocal
654 c_nl => spotential(kkind)%sgp_potential%c_nonlocal
655 ppnl_radius = spotential(kkind)%sgp_potential%ppnl_radius
656 ALLOCATE (radp(nnl))
657 radp(:) = ppnl_radius
658 cprj => spotential(kkind)%sgp_potential%cprj_ppnl
659 hprj => spotential(kkind)%sgp_potential%vprj_ppnl
660 nppnl = SIZE(cprj, 2)
661 END IF
662 ELSE
663 cycle
664 END IF
665
666 IF (my_ref) THEN
667 raf(:) = pbc(particle_set(iatom)%r(:) - rf, cell)
668 ELSE
669 raf(:) = -rac
670 END IF
671 ! Keep both moment centers on the image used by the projector neighbor list.
672 rcf(:) = raf + rac
673
674 dac = norm2(rac)
675 clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
676 clist%catom = katom
677 clist%cell = cell_c
678 clist%rac = rac
679 ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
680 clist%achint(nsgfa, nppnl, maxder))
681 clist%acint = 0._dp
682 clist%achint = 0._dp
683 clist%nsgf_cnt = 0
684 NULLIFY (clist%sgf_list)
685
686 IF (PRESENT(pseudoatom)) THEN
687 IF (pseudoatom /= katom) THEN
688 clist%maxac = maxval(abs(clist%acint(:, :, 1)))
689 clist%maxach = maxval(abs(clist%achint(:, :, 1)))
690 IF (.NOT. dogth .AND. .NOT. do_upf) DEALLOCATE (radp)
691 cycle
692 END IF
693 END IF
694
695 DO iset = 1, nseta
696 ncoa = npgfa(iset)*ncoset(la_max(iset))
697 sgfa = first_sgfa(1, iset)
698 IF (dogth) THEN
699 ! GTH potential
700 prjc = 1
701 work = 0._dp
702 DO l = 0, lppnl
703 nprjc = nprj_ppnl(l)*nco(l)
704 IF (nprjc == 0) cycle
705 rprjc(1) = ppnl_radius
706 IF (set_radius_a(iset) + rprjc(1) < dac) cycle
707 lc_max = l + 2*(nprj_ppnl(l) - 1)
708 lc_min = l
709 zetc(1) = alpha_ppnl(l)
710 ncoc = ncoset(lc_max)
711
712 ! *** Calculate the primitive overlap integrals ***
713 IF (my_momentmode) THEN
714 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
715 lc_max, 1, zetc, rprjc, nder, raf, rcf, ai_work, overlap=sab)
716 ! reduce ai_work to sab
717 na = ncoa
718 np = ncoc
719 DO i = 1, maxder - 1
720 first_col = i*ldsab
721 sab(1:na, first_col + 1:first_col + np) = ai_work(1:na, 1:np, i)
722 END DO
723 ELSE
724 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
725 lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .true., ai_work, ldai)
726 END IF
727
728 ! *** Transformation step projector functions (cartesian->spherical) ***
729 na = ncoa
730 nb = nprjc
731 np = ncoc
732 DO i = 1, maxder
733 first_col = (i - 1)*ldsab
734 work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
735 matmul(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
736 END DO
737 prjc = prjc + nprjc
738 END DO
739
740 na = nsgf_seta(iset)
741 nb = nppnl
742 np = ncoa
743 DO i = 1, maxder
744 first_col = (i - 1)*ldsab + 1
745
746 ! *** Contraction step (basis functions) ***
747 clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
748 matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))
749
750 ! *** Multiply with interaction matrix(h) ***
751 clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
752 matmul(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))
753 END DO
754 ELSE IF (do_upf) THEN
755 na = nsgf_seta(iset)
756 nb = nppnl
757 np = ncoa
758 ALLOCATE (upf_int(np, nb, maxder))
759 IF (my_momentmode) THEN
760 CALL upf_projector_overlap(upf, la_max(iset), zeta(:npgfa(iset), iset), &
761 rac, nder, upf_int, moment_origin=raf)
762 ELSE
763 CALL upf_projector_overlap(upf, la_max(iset), zeta(:npgfa(iset), iset), &
764 rac, nder, upf_int)
765 END IF
766 DO i = 1, maxder
767 clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
768 matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), upf_int(:, :, i))
769 clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
770 matmul(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), upf%d)
771 END DO
772 DEALLOCATE (upf_int)
773 ELSE
774 ! SGP potential
775 ! *** Calculate the primitive overlap integrals ***
776 IF (my_momentmode) THEN
777 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
778 lppnl, nnl, a_nl, radp, nder, raf, rcf, ai_work, overlap=sab)
779 ! reduce ai_work to sab
780 na = ncoa
781 DO i = 1, maxder - 1
782 first_col = i*ldsab + 1
783 sab(1:na, first_col:first_col + nprjc - 1) = ai_work(1:na, 1:nprjc, i)
784 END DO
785 ELSE
786 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
787 lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .true., ai_work, ldai)
788 END IF
789
790 na = nsgf_seta(iset)
791 nb = nppnl
792 np = ncoa
793 DO i = 1, maxder
794 first_col = (i - 1)*ldsab + 1
795 ! *** Transformation step projector functions (cartesian->spherical) ***
796 work(1:np, 1:nb) = matmul(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb))
797
798 ! *** Contraction step (basis functions) ***
799 clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
800 matmul(transpose(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb))
801
802 ! *** Multiply with interaction matrix(h) ***
803 ncoc = sgfa + nsgf_seta(iset) - 1
804 DO j = 1, nppnl
805 clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
806 END DO
807 END DO
808 END IF
809 END DO
810 clist%maxac = maxval(abs(clist%acint(:, :, 1)))
811 clist%maxach = maxval(abs(clist%achint(:, :, 1)))
812 IF (.NOT. dogth .AND. .NOT. do_upf) DEALLOCATE (radp)
813
814 END DO
815
816 DEALLOCATE (sab, ai_work, work)
817
818!$OMP END PARALLEL
819
820 DEALLOCATE (basis_set, gpotential, spotential)
821
822 CALL timestop(handle)
823
824 END SUBROUTINE build_sap_ints
825
826END MODULE sap_kind_types
cell_types::pbc
Definition cell_types.F:105
sap_kind_types::get_alist
Definition sap_kind_types.F:113
sap_kind_types::release_sap_int
Definition sap_kind_types.F:110
sap_kind_types::sap_sort
Definition sap_kind_types.F:116
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
ai_moments::moment
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab, overlap)
...
Definition ai_moments.F:289
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition ai_overlap.F:75
atom
Definition atom.F:9
basis_set_types
Definition basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_kind_types::get_qs_kind_set
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, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
sap_kind_types
General overlap type integrals containers.
Definition sap_kind_types.F:12
sap_kind_types::alist_post_align_blk
subroutine, public alist_post_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
...
Definition sap_kind_types.F:333
sap_kind_types::build_sap_ints
subroutine, public build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder, moment_mode, refpoint, particle_set, cell, pseudoatom)
Calculate overlap and optionally momenta <a|x^n|p> between GTOs and nl. pseudo potential projectors a...
Definition sap_kind_types.F:425
sap_kind_types::alist_pre_align_blk
subroutine, public alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
...
Definition sap_kind_types.F:292
upf_projector_integrals
Half-numerical overlaps with the original scalar UPF radial projectors. Cartesian polynomials and ang...
Definition upf_projector_integrals.F:15
upf_projector_integrals::upf_projector_init
subroutine, public upf_projector_init(pot, degree, projectors)
Prepare immutable radial and angular data outside parallel regions.
Definition upf_projector_integrals.F:64
upf_projector_integrals::upf_projector_overlap
subroutine, public upf_projector_overlap(projectors, la, zeta, rac, nder, integrals, moment_origin)
Evaluate <Cartesian Gaussian|beta_i Y_lm> and AO-center derivatives.
Definition upf_projector_integrals.F:184
util
All kind of helpful little routines.
Definition util.F:14
util::locate
pure integer function, public locate(array, x)
Purpose: Given an array array(1:n), and given a value x, a value x_index is returned which is the ind...
Definition util.F:61
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
external_potential_types::gth_potential_p_type
Definition external_potential_types.F:238
external_potential_types::gth_potential_type
Definition external_potential_types.F:121
external_potential_types::sgp_potential_p_type
Definition external_potential_types.F:246
external_potential_types::sgp_potential_type
Definition external_potential_types.F:179
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
sap_kind_types::alist_complex_type
Definition sap_kind_types.F:90
sap_kind_types::alist_type
Definition sap_kind_types.F:61
sap_kind_types::clist_complex_type
Definition sap_kind_types.F:75
sap_kind_types::clist_type
Definition sap_kind_types.F:46
sap_kind_types::sap_int_complex_type
Definition sap_kind_types.F:96
sap_kind_types::sap_int_type
Definition sap_kind_types.F:67
upf_projector_integrals::upf_projector_type
Definition upf_projector_integrals.F:37