Line data Source code
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 : ! **************************************************************************************************
12 : MODULE 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 : #:for type_suffix, integral_type in [('', 'REAL'), ('_complex', 'COMPLEX')]
47 : TYPE clist${type_suffix}$_type
48 : INTEGER :: catom = -1
49 : INTEGER :: nsgf_cnt = -1
50 : INTEGER, DIMENSION(:), POINTER :: sgf_list => NULL()
51 : INTEGER, DIMENSION(3) :: cell = -1
52 : LOGICAL :: sgf_soft_only = .FALSE.
53 : REAL(KIND=dp) :: maxac = 0.0_dp
54 : REAL(KIND=dp) :: maxach = 0.0_dp
55 : REAL(KIND=dp), DIMENSION(3) :: rac = 0.0_dp
56 : ${integral_type}$ (KIND=dp), DIMENSION(:, :, :), POINTER :: acint => NULL()
57 : ${integral_type}$ (KIND=dp), DIMENSION(:, :, :), POINTER :: achint => NULL()
58 : ${integral_type}$ (KIND=dp), DIMENSION(:, :, :), POINTER :: alint => NULL()
59 : ${integral_type}$ (KIND=dp), DIMENSION(:, :, :), POINTER :: alkint => NULL()
60 : END TYPE clist${type_suffix}$_type
61 :
62 : TYPE alist${type_suffix}$_type
63 : INTEGER :: aatom = -1
64 : INTEGER :: nclist = -1
65 : TYPE(clist${type_suffix}$_type), DIMENSION(:), POINTER :: clist => NULL()
66 : END TYPE alist${type_suffix}$_type
67 :
68 : TYPE sap_int${type_suffix}$_type
69 : INTEGER :: a_kind = -1
70 : INTEGER :: p_kind = -1
71 : INTEGER :: nalist = -1
72 : TYPE(alist${type_suffix}$_type), DIMENSION(:), POINTER :: alist => NULL()
73 : INTEGER, DIMENSION(:), POINTER :: asort => NULL()
74 : INTEGER, DIMENSION(:), POINTER :: aindex => NULL()
75 : END TYPE sap_int${type_suffix}$_type
76 : #:endfor
77 :
78 : PUBLIC :: sap_int_type, clist_type, alist_type, &
79 : sap_int_complex_type, clist_complex_type, alist_complex_type, &
80 : release_sap_int, get_alist, alist_pre_align_blk, &
81 : alist_post_align_blk, sap_sort, build_sap_ints
82 :
83 : INTERFACE release_sap_int
84 : MODULE PROCEDURE release_sap_int_real, release_sap_int_complex
85 : END INTERFACE release_sap_int
86 : INTERFACE get_alist
87 : MODULE PROCEDURE get_alist_real, get_alist_complex
88 : END INTERFACE get_alist
89 : INTERFACE sap_sort
90 : MODULE PROCEDURE sap_sort_real, sap_sort_complex
91 : END INTERFACE sap_sort
92 :
93 : CONTAINS
94 :
95 : !==========================================================================================================
96 :
97 : #:for type_suffix, suffix in [('', 'real'), ('_complex', 'complex')]
98 : ! **************************************************************************************************
99 : !> \brief ...
100 : !> \param sap_int ...
101 : ! **************************************************************************************************
102 22034 : SUBROUTINE release_sap_int_${suffix}$ (sap_int)
103 :
104 : TYPE(sap_int${type_suffix}$_type), DIMENSION(:), POINTER :: sap_int
105 :
106 : INTEGER :: i, j, k
107 : TYPE(clist${type_suffix}$_type), POINTER :: clist
108 :
109 22034 : CPASSERT(ASSOCIATED(sap_int))
110 :
111 107428 : DO i = 1, SIZE(sap_int)
112 85394 : IF (ASSOCIATED(sap_int(i)%alist)) THEN
113 159420 : DO j = 1, SIZE(sap_int(i)%alist)
114 159420 : IF (ASSOCIATED(sap_int(i)%alist(j)%clist)) THEN
115 1710185 : DO k = 1, SIZE(sap_int(i)%alist(j)%clist)
116 1602745 : clist => sap_int(i)%alist(j)%clist(k)
117 1602745 : IF (ASSOCIATED(clist%acint)) THEN
118 1586573 : DEALLOCATE (clist%acint)
119 : END IF
120 1602745 : IF (ASSOCIATED(clist%sgf_list)) THEN
121 137350 : DEALLOCATE (clist%sgf_list)
122 : END IF
123 1602745 : IF (ASSOCIATED(clist%achint)) THEN
124 1157622 : DEALLOCATE (clist%achint)
125 : END IF
126 1602745 : IF (ASSOCIATED(clist%alint)) THEN
127 1138244 : DEALLOCATE (clist%alint)
128 : END IF
129 1710185 : IF (ASSOCIATED(clist%alkint)) THEN
130 1138244 : DEALLOCATE (clist%alkint)
131 : END IF
132 : END DO
133 107440 : DEALLOCATE (sap_int(i)%alist(j)%clist)
134 : END IF
135 : END DO
136 51897 : DEALLOCATE (sap_int(i)%alist)
137 : END IF
138 85394 : IF (ASSOCIATED(sap_int(i)%asort)) THEN
139 51114 : DEALLOCATE (sap_int(i)%asort)
140 : END IF
141 107428 : IF (ASSOCIATED(sap_int(i)%aindex)) THEN
142 51114 : DEALLOCATE (sap_int(i)%aindex)
143 : END IF
144 : END DO
145 :
146 22034 : DEALLOCATE (sap_int)
147 :
148 22034 : END SUBROUTINE release_sap_int_${suffix}$
149 : #:endfor
150 :
151 : #:for type_suffix, suffix in [('', 'real'), ('_complex', 'complex')]
152 : ! **************************************************************************************************
153 : !> \brief ...
154 : !> \param sap_int ...
155 : !> \param alist ...
156 : !> \param atom ...
157 : ! **************************************************************************************************
158 10996970 : SUBROUTINE get_alist_${suffix}$ (sap_int, alist, atom)
159 :
160 : TYPE(sap_int${type_suffix}$_type), INTENT(IN) :: sap_int
161 : TYPE(alist${type_suffix}$_type), INTENT(OUT), POINTER :: alist
162 : INTEGER, INTENT(IN) :: atom
163 :
164 : INTEGER :: i
165 :
166 10996970 : NULLIFY (alist)
167 10996970 : i = locate(sap_int%asort, atom)
168 10996970 : IF (i > 0 .AND. i <= SIZE(sap_int%alist)) THEN
169 10995909 : i = sap_int%aindex(i)
170 10995909 : alist => sap_int%alist(i)
171 1061 : ELSE IF (i == 0) THEN
172 : NULLIFY (alist)
173 : ELSE
174 0 : CPABORT("Invalid result of locate(sap_int%asort, atom)")
175 : END IF
176 :
177 10996970 : END SUBROUTINE get_alist_${suffix}$
178 : #:endfor
179 :
180 : ! **************************************************************************************************
181 : !> \brief ...
182 : !> \param blk_in ...
183 : !> \param ldin ...
184 : !> \param blk_out ...
185 : !> \param ldout ...
186 : !> \param ilist ...
187 : !> \param in ...
188 : !> \param jlist ...
189 : !> \param jn ...
190 : ! **************************************************************************************************
191 8877084 : SUBROUTINE alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
192 : INTEGER, INTENT(IN) :: in, ilist(*), ldout
193 : REAL(dp), INTENT(INOUT) :: blk_out(ldout, *)
194 : INTEGER, INTENT(IN) :: ldin
195 : REAL(dp), INTENT(IN) :: blk_in(ldin, *)
196 : INTEGER, INTENT(IN) :: jlist(*), jn
197 :
198 : INTEGER :: i, i0, i1, i2, i3, inn, inn1, j, j0
199 :
200 8877084 : inn = MOD(in, 4)
201 8877084 : inn1 = inn + 1
202 63831898 : DO j = 1, jn
203 54954814 : j0 = jlist(j)
204 88565605 : DO i = 1, inn
205 33610791 : i0 = ilist(i)
206 88565605 : blk_out(i, j) = blk_in(i0, j0)
207 : END DO
208 63831898 : DO i = inn1, in, 4
209 97718789 : i0 = ilist(i)
210 97718789 : i1 = ilist(i + 1)
211 97718789 : i2 = ilist(i + 2)
212 97718789 : i3 = ilist(i + 3)
213 97718789 : blk_out(i, j) = blk_in(i0, j0)
214 97718789 : blk_out(i + 1, j) = blk_in(i1, j0)
215 97718789 : blk_out(i + 2, j) = blk_in(i2, j0)
216 104053394 : blk_out(i + 3, j) = blk_in(i3, j0)
217 : END DO
218 : END DO
219 8877084 : END SUBROUTINE alist_pre_align_blk
220 :
221 : ! **************************************************************************************************
222 : !> \brief ...
223 : !> \param blk_in ...
224 : !> \param ldin ...
225 : !> \param blk_out ...
226 : !> \param ldout ...
227 : !> \param ilist ...
228 : !> \param in ...
229 : !> \param jlist ...
230 : !> \param jn ...
231 : ! **************************************************************************************************
232 8621703 : SUBROUTINE alist_post_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
233 : INTEGER, INTENT(IN) :: in, ilist(*), ldout
234 : REAL(dp), INTENT(INOUT) :: blk_out(ldout, *)
235 : INTEGER, INTENT(IN) :: ldin
236 : REAL(dp), INTENT(IN) :: blk_in(ldin, *)
237 : INTEGER, INTENT(IN) :: jlist(*), jn
238 :
239 : INTEGER :: i, i0, i1, i2, i3, inn, inn1, j, j0
240 :
241 8621703 : inn = MOD(in, 4)
242 8621703 : inn1 = inn + 1
243 61733142 : DO j = 1, jn
244 53111439 : j0 = jlist(j)
245 86465128 : DO i = 1, inn
246 33353689 : i0 = ilist(i)
247 86465128 : blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j)
248 : END DO
249 61733142 : DO i = inn1, in, 4
250 93349830 : i0 = ilist(i)
251 93349830 : i1 = ilist(i + 1)
252 93349830 : i2 = ilist(i + 2)
253 93349830 : i3 = ilist(i + 3)
254 93349830 : blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j)
255 93349830 : blk_out(i1, j0) = blk_out(i1, j0) + blk_in(i + 1, j)
256 93349830 : blk_out(i2, j0) = blk_out(i2, j0) + blk_in(i + 2, j)
257 99263226 : blk_out(i3, j0) = blk_out(i3, j0) + blk_in(i + 3, j)
258 : END DO
259 : END DO
260 8621703 : END SUBROUTINE alist_post_align_blk
261 :
262 : #:for type_suffix, suffix in [('', 'real'), ('_complex', 'complex')]
263 : ! **************************************************************************************************
264 : !> \brief ...
265 : !> \param sap_int ...
266 : ! **************************************************************************************************
267 21660 : SUBROUTINE sap_sort_${suffix}$ (sap_int)
268 : TYPE(sap_int${type_suffix}$_type), DIMENSION(:), POINTER :: sap_int
269 :
270 : INTEGER :: iac, na
271 :
272 : ! *** Set up a sorting index
273 :
274 21660 : !$OMP PARALLEL DEFAULT(NONE) SHARED(sap_int) PRIVATE(iac,na)
275 : !$OMP DO
276 : DO iac = 1, SIZE(sap_int)
277 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
278 : na = SIZE(sap_int(iac)%alist)
279 : ALLOCATE (sap_int(iac)%asort(na), sap_int(iac)%aindex(na))
280 : sap_int(iac)%asort(1:na) = sap_int(iac)%alist(1:na)%aatom
281 : CALL sort(sap_int(iac)%asort, na, sap_int(iac)%aindex)
282 : END DO
283 : !$OMP END PARALLEL
284 :
285 21660 : END SUBROUTINE sap_sort_${suffix}$
286 : #:endfor
287 :
288 : !==========================================================================================================
289 :
290 : ! **************************************************************************************************
291 : !> \brief Calculate overlap and optionally momenta <a|x^n|p> between GTOs and nl. pseudo potential projectors
292 : !> adapted from core_ppnl.F::build_core_ppnl
293 : !> \param sap_int allocated in parent routine (nkind*nkind)
294 : !> \param sap_ppnl ...
295 : !> \param qs_kind_set ...
296 : !> \param nder Either number of derivatives or order of moments
297 : !> \param moment_mode if present and true, moments are calculated instead of derivatives
298 : !> \param refpoint optionally the reference point for moment calculation
299 : !> \param particle_set needed if refpoint is present
300 : !> \param cell needed if refpoint is present
301 : !> \param pseudoatom If we want to constrain the calculations to katom == pseudoatom
302 : ! **************************************************************************************************
303 134 : SUBROUTINE build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder, moment_mode, refpoint, particle_set, cell, pseudoatom)
304 : TYPE(sap_int_type), DIMENSION(:), INTENT(INOUT), &
305 : POINTER :: sap_int
306 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
307 : INTENT(IN), POINTER :: sap_ppnl
308 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
309 : POINTER :: qs_kind_set
310 : INTEGER, INTENT(IN) :: nder
311 : LOGICAL, INTENT(IN), OPTIONAL :: moment_mode
312 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
313 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
314 : OPTIONAL, POINTER :: particle_set
315 : TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
316 : INTEGER, OPTIONAL :: pseudoatom
317 :
318 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_sap_ints'
319 :
320 : INTEGER :: first_col, handle, i, iac, iatom, ikind, ilist, iset, j, jneighbor, katom, kkind, &
321 : l, lc_max, lc_min, ldai, ldai_nl, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, &
322 : maxppnl, maxsgf, na, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, &
323 : nseta, nsgfa, prjc, sgfa, slot
324 : INTEGER, DIMENSION(3) :: cell_c
325 134 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
326 134 : nsgf_seta
327 134 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
328 : LOGICAL :: do_upf, dogth, my_momentmode, my_ref
329 : LOGICAL, DIMENSION(0:9) :: is_nonlocal
330 : REAL(KIND=dp) :: dac, ppnl_radius
331 134 : REAL(KIND=dp), ALLOCATABLE :: upf_int(:, :, :)
332 134 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp
333 134 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
334 134 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work
335 : REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
336 : REAL(KIND=dp), DIMENSION(3) :: ra, rac, raf, rc, rcf, rf
337 134 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
338 134 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_nl, rpgfa, sphi_a, vprj_ppnl, &
339 134 : zeta
340 134 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nl
341 : TYPE(clist_type), POINTER :: clist
342 134 : TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
343 : TYPE(gth_potential_type), POINTER :: gth_potential
344 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
345 134 : DIMENSION(:) :: basis_set
346 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
347 134 : TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
348 : TYPE(sgp_potential_type), POINTER :: sgp_potential
349 : TYPE(upf_projector_type), POINTER :: upf
350 :
351 134 : CALL timeset(routineN, handle)
352 :
353 134 : nkind = SIZE(qs_kind_set)
354 134 : maxder = ncoset(nder)
355 :
356 : ! determine whether moments or derivatives should be calculated
357 134 : my_momentmode = .FALSE.
358 134 : IF (PRESENT(moment_mode)) THEN
359 134 : my_momentmode = moment_mode
360 : END IF
361 :
362 134 : my_ref = .FALSE.
363 134 : IF (PRESENT(refpoint)) THEN
364 : ! A reference point also requires the cell and particle coordinates.
365 64 : CPASSERT((PRESENT(cell) .AND. PRESENT(particle_set)))
366 64 : rf = refpoint
367 64 : my_ref = .TRUE.
368 : END IF
369 :
370 : CALL get_qs_kind_set(qs_kind_set, &
371 : maxco=maxco, &
372 : maxlgto=maxlgto, &
373 : maxsgf=maxsgf, &
374 : maxlppnl=maxlppnl, &
375 134 : maxppnl=maxppnl)
376 :
377 134 : maxl = MAX(maxlgto, maxlppnl)
378 134 : CALL init_orbital_pointers(maxl + nder + 1)
379 :
380 134 : ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
381 134 : IF (.NOT. my_momentmode) THEN
382 0 : ldai = ncoset(maxl + nder + 1)
383 : ELSE
384 134 : ldai = maxco
385 : END IF
386 :
387 : !set up direct access to basis and potential
388 134 : ldai_nl = 0
389 134 : NULLIFY (gpotential, spotential)
390 1570 : ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
391 434 : DO ikind = 1, nkind
392 300 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
393 300 : IF (ASSOCIATED(orb_basis_set)) THEN
394 300 : basis_set(ikind)%gto_basis_set => orb_basis_set
395 : ELSE
396 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
397 : END IF
398 300 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
399 300 : NULLIFY (gpotential(ikind)%gth_potential)
400 300 : NULLIFY (spotential(ikind)%sgp_potential)
401 434 : IF (ASSOCIATED(gth_potential)) THEN
402 300 : gpotential(ikind)%gth_potential => gth_potential
403 300 : ldai_nl = MAX(ldai_nl, ncoset(gth_potential%lprj_ppnl_max))
404 0 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
405 0 : spotential(ikind)%sgp_potential => sgp_potential
406 0 : IF (ALLOCATED(sgp_potential%upf_data)) THEN
407 0 : IF (sgp_potential%upf_data%number_of_proj > 0) THEN
408 0 : IF (sgp_potential%upf_projectors%degree < maxlgto + nder) THEN
409 0 : CALL upf_projector_init(sgp_potential%upf_data, maxlgto + nder, sgp_potential%upf_projectors)
410 : END IF
411 : END IF
412 : END IF
413 0 : ldai_nl = MAX(ldai_nl, sgp_potential%n_nonlocal*ncoset(sgp_potential%lmax))
414 : END IF
415 : END DO
416 :
417 : !allocate sap int
418 6288 : DO slot = 1, sap_ppnl(1)%nl_size
419 :
420 6154 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
421 6154 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
422 6154 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
423 6154 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
424 6154 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
425 6154 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
426 6154 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
427 :
428 6154 : iac = ikind + nkind*(kkind - 1)
429 6154 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
430 6154 : IF (.NOT. ASSOCIATED(gpotential(kkind)%gth_potential) .AND. &
431 : .NOT. ASSOCIATED(spotential(kkind)%sgp_potential)) CYCLE
432 6154 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
433 476 : sap_int(iac)%a_kind = ikind
434 476 : sap_int(iac)%p_kind = kkind
435 476 : sap_int(iac)%nalist = nlist
436 3534 : ALLOCATE (sap_int(iac)%alist(nlist))
437 2582 : DO i = 1, nlist
438 2106 : NULLIFY (sap_int(iac)%alist(i)%clist)
439 2106 : sap_int(iac)%alist(i)%aatom = 0
440 2582 : sap_int(iac)%alist(i)%nclist = 0
441 : END DO
442 : END IF
443 6288 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
444 2106 : sap_int(iac)%alist(ilist)%aatom = iatom
445 2106 : sap_int(iac)%alist(ilist)%nclist = nneighbor
446 27214 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
447 8260 : DO i = 1, nneighbor
448 8260 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
449 : END DO
450 : END IF
451 : END DO
452 :
453 : !calculate the overlap integrals <a|pp>
454 : !$OMP PARALLEL &
455 : !$OMP DEFAULT (NONE) &
456 : !$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, my_momentmode, ldai_nl, &
457 : !$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco, my_ref, cell, particle_set, rf, pseudoatom) &
458 : !$OMP PRIVATE (ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
459 : !$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
460 : !$OMP slot, sphi_a, zeta, cprj, hprj, lppnl, nppnl, nprj_ppnl, &
461 : !$OMP do_upf, upf, upf_int, &
462 : !$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, &
463 : !$OMP ncoc, rpgfa, first_col, vprj_ppnl, i, j, l, dogth, &
464 : !$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl, &
465 134 : !$OMP na, nb, np, nnl, is_nonlocal, a_nl, h_nl, c_nl, radp, raf, rcf, ra, rc)
466 :
467 : ! allocate work storage
468 : ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
469 : sab = 0.0_dp
470 : IF (.NOT. my_momentmode) THEN
471 : ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
472 : ELSE
473 : ALLOCATE (ai_work(ldai, ldai_nl, ncoset(nder + 1)))
474 : END IF
475 : ai_work = 0.0_dp
476 :
477 : !$OMP DO SCHEDULE(GUIDED)
478 : ! loop over neighbourlist
479 : DO slot = 1, sap_ppnl(1)%nl_size
480 : ikind = sap_ppnl(1)%nlist_task(slot)%ikind
481 : kkind = sap_ppnl(1)%nlist_task(slot)%jkind
482 : iatom = sap_ppnl(1)%nlist_task(slot)%iatom
483 : katom = sap_ppnl(1)%nlist_task(slot)%jatom
484 : nlist = sap_ppnl(1)%nlist_task(slot)%nlist
485 : ilist = sap_ppnl(1)%nlist_task(slot)%ilist
486 : nneighbor = sap_ppnl(1)%nlist_task(slot)%nnode
487 : jneighbor = sap_ppnl(1)%nlist_task(slot)%inode
488 : cell_c(:) = sap_ppnl(1)%nlist_task(slot)%cell(:)
489 : rac(1:3) = sap_ppnl(1)%nlist_task(slot)%r(1:3)
490 :
491 : iac = ikind + nkind*(kkind - 1)
492 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
493 : ! get definition of basis set
494 : first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
495 : la_max => basis_set(ikind)%gto_basis_set%lmax
496 : la_min => basis_set(ikind)%gto_basis_set%lmin
497 : npgfa => basis_set(ikind)%gto_basis_set%npgf
498 : nseta = basis_set(ikind)%gto_basis_set%nset
499 : nsgfa = basis_set(ikind)%gto_basis_set%nsgf
500 : nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
501 : rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
502 : set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
503 : sphi_a => basis_set(ikind)%gto_basis_set%sphi
504 : zeta => basis_set(ikind)%gto_basis_set%zet
505 : ! get definition of PP projectors
506 : do_upf = .FALSE.
507 : IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN
508 : ! GTH potential
509 : dogth = .TRUE.
510 : alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
511 : cprj => gpotential(kkind)%gth_potential%cprj
512 : lppnl = gpotential(kkind)%gth_potential%lppnl
513 : nppnl = gpotential(kkind)%gth_potential%nppnl
514 : nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
515 : ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
516 : vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
517 : ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN
518 : ! SGP potential
519 : dogth = .FALSE.
520 : nprjc = spotential(kkind)%sgp_potential%nppnl
521 : IF (nprjc == 0) CYCLE
522 : do_upf = ALLOCATED(spotential(kkind)%sgp_potential%upf_data)
523 : IF (do_upf) THEN
524 : upf => spotential(kkind)%sgp_potential%upf_projectors
525 : nppnl = upf%nproj
526 : ELSE
527 : nnl = spotential(kkind)%sgp_potential%n_nonlocal
528 : lppnl = spotential(kkind)%sgp_potential%lmax
529 : is_nonlocal = .FALSE.
530 : is_nonlocal(0:lppnl) = spotential(kkind)%sgp_potential%is_nonlocal(0:lppnl)
531 : a_nl => spotential(kkind)%sgp_potential%a_nonlocal
532 : h_nl => spotential(kkind)%sgp_potential%h_nonlocal
533 : c_nl => spotential(kkind)%sgp_potential%c_nonlocal
534 : ppnl_radius = spotential(kkind)%sgp_potential%ppnl_radius
535 : ALLOCATE (radp(nnl))
536 : radp(:) = ppnl_radius
537 : cprj => spotential(kkind)%sgp_potential%cprj_ppnl
538 : hprj => spotential(kkind)%sgp_potential%vprj_ppnl
539 : nppnl = SIZE(cprj, 2)
540 : END IF
541 : ELSE
542 : CYCLE
543 : END IF
544 :
545 : IF (my_ref) THEN
546 : ra(:) = pbc(particle_set(iatom)%r(:) - rf, cell) + rf
547 : rc(:) = pbc(particle_set(katom)%r(:) - rf, cell) + rf
548 : raf(:) = ra(:) - rf(:)
549 : rcf(:) = rc(:) - rf(:)
550 : ELSE
551 : raf(:) = -rac
552 : rcf(:) = [0._dp, 0._dp, 0._dp]
553 : END IF
554 :
555 : dac = NORM2(rac)
556 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
557 : clist%catom = katom
558 : clist%cell = cell_c
559 : clist%rac = rac
560 : ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
561 : clist%achint(nsgfa, nppnl, maxder))
562 : clist%acint = 0._dp
563 : clist%achint = 0._dp
564 : clist%nsgf_cnt = 0
565 : NULLIFY (clist%sgf_list)
566 :
567 : IF (PRESENT(pseudoatom)) THEN
568 : IF (pseudoatom /= katom) THEN
569 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
570 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
571 : IF (.NOT. dogth .AND. .NOT. do_upf) DEALLOCATE (radp)
572 : CYCLE
573 : END IF
574 : END IF
575 :
576 : DO iset = 1, nseta
577 : ncoa = npgfa(iset)*ncoset(la_max(iset))
578 : sgfa = first_sgfa(1, iset)
579 : IF (dogth) THEN
580 : ! GTH potential
581 : prjc = 1
582 : work = 0._dp
583 : DO l = 0, lppnl
584 : nprjc = nprj_ppnl(l)*nco(l)
585 : IF (nprjc == 0) CYCLE
586 : rprjc(1) = ppnl_radius
587 : IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
588 : lc_max = l + 2*(nprj_ppnl(l) - 1)
589 : lc_min = l
590 : zetc(1) = alpha_ppnl(l)
591 : ncoc = ncoset(lc_max)
592 :
593 : ! *** Calculate the primitive overlap integrals ***
594 : IF (my_momentmode) THEN
595 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
596 : lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, 0, .FALSE., ai_work, ldai)
597 : ELSE
598 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
599 : lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .TRUE., ai_work, ldai)
600 : END IF
601 : IF (my_momentmode .AND. nder >= 1) THEN
602 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
603 : lc_max, 1, zetc, rprjc, nder, raf, rcf, ai_work)
604 : ! reduce ai_work to sab
605 : na = ncoa
606 : np = ncoc
607 : DO i = 1, maxder - 1
608 : first_col = i*ldsab
609 : sab(1:na, first_col + 1:first_col + np) = ai_work(1:na, 1:np, i)
610 : END DO
611 : END IF
612 :
613 : ! *** Transformation step projector functions (cartesian->spherical) ***
614 : na = ncoa
615 : nb = nprjc
616 : np = ncoc
617 : DO i = 1, maxder
618 : first_col = (i - 1)*ldsab
619 : work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
620 : MATMUL(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
621 : END DO
622 : prjc = prjc + nprjc
623 : END DO
624 :
625 : na = nsgf_seta(iset)
626 : nb = nppnl
627 : np = ncoa
628 : DO i = 1, maxder
629 : first_col = (i - 1)*ldsab + 1
630 :
631 : ! *** Contraction step (basis functions) ***
632 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
633 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))
634 :
635 : ! *** Multiply with interaction matrix(h) ***
636 : clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
637 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))
638 : END DO
639 : ELSE IF (do_upf) THEN
640 : na = nsgf_seta(iset)
641 : nb = nppnl
642 : np = ncoa
643 : ALLOCATE (upf_int(np, nb, maxder))
644 : IF (my_momentmode) THEN
645 : CALL upf_projector_overlap(upf, la_max(iset), zeta(:npgfa(iset), iset), &
646 : rac, nder, upf_int, moment_origin=raf)
647 : ELSE
648 : CALL upf_projector_overlap(upf, la_max(iset), zeta(:npgfa(iset), iset), &
649 : rac, nder, upf_int)
650 : END IF
651 : DO i = 1, maxder
652 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
653 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), upf_int(:, :, i))
654 : clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
655 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), upf%d)
656 : END DO
657 : DEALLOCATE (upf_int)
658 : ELSE
659 : ! SGP potential
660 : ! *** Calculate the primitive overlap integrals ***
661 : IF (my_momentmode) THEN
662 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
663 : lppnl, 0, nnl, radp, a_nl, rac, dac, sab, 0, .FALSE., ai_work, ldai)
664 : ELSE
665 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
666 : lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .TRUE., ai_work, ldai)
667 : END IF
668 : IF (my_momentmode .AND. nder >= 1) THEN
669 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
670 : lppnl, nnl, a_nl, radp, nder, raf, rcf, ai_work)
671 : ! reduce ai_work to sab
672 : na = ncoa
673 : DO i = 1, maxder - 1
674 : first_col = i*ldsab
675 : sab(1:na, first_col:first_col + nprjc - 1) = ai_work(1:na, 1:nprjc, i)
676 : END DO
677 : END IF
678 :
679 : na = nsgf_seta(iset)
680 : nb = nppnl
681 : np = ncoa
682 : DO i = 1, maxder
683 : first_col = (i - 1)*ldsab + 1
684 : ! *** Transformation step projector functions (cartesian->spherical) ***
685 : work(1:np, 1:nb) = MATMUL(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb))
686 :
687 : ! *** Contraction step (basis functions) ***
688 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
689 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb))
690 :
691 : ! *** Multiply with interaction matrix(h) ***
692 : ncoc = sgfa + nsgf_seta(iset) - 1
693 : DO j = 1, nppnl
694 : clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
695 : END DO
696 : END DO
697 : END IF
698 : END DO
699 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
700 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
701 : IF (.NOT. dogth .AND. .NOT. do_upf) DEALLOCATE (radp)
702 :
703 : END DO
704 :
705 : DEALLOCATE (sab, ai_work, work)
706 :
707 : !$OMP END PARALLEL
708 :
709 134 : DEALLOCATE (basis_set, gpotential, spotential)
710 :
711 134 : CALL timestop(handle)
712 :
713 402 : END SUBROUTINE build_sap_ints
714 :
715 0 : END MODULE sap_kind_types
|
