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 22082 : 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 22082 : CPASSERT(ASSOCIATED(sap_int))
110 :
111 107608 : DO i = 1, SIZE(sap_int)
112 85526 : IF (ASSOCIATED(sap_int(i)%alist)) THEN
113 159938 : DO j = 1, SIZE(sap_int(i)%alist)
114 159938 : IF (ASSOCIATED(sap_int(i)%alist(j)%clist)) THEN
115 1714131 : DO k = 1, SIZE(sap_int(i)%alist(j)%clist)
116 1606301 : clist => sap_int(i)%alist(j)%clist(k)
117 1606301 : IF (ASSOCIATED(clist%acint)) THEN
118 1590129 : DEALLOCATE (clist%acint)
119 : END IF
120 1606301 : IF (ASSOCIATED(clist%sgf_list)) THEN
121 137438 : DEALLOCATE (clist%sgf_list)
122 : END IF
123 1606301 : IF (ASSOCIATED(clist%achint)) THEN
124 1161114 : DEALLOCATE (clist%achint)
125 : END IF
126 1606301 : IF (ASSOCIATED(clist%alint)) THEN
127 1141712 : DEALLOCATE (clist%alint)
128 : END IF
129 1714131 : IF (ASSOCIATED(clist%alkint)) THEN
130 1141712 : DEALLOCATE (clist%alkint)
131 : END IF
132 : END DO
133 107830 : DEALLOCATE (sap_int(i)%alist(j)%clist)
134 : END IF
135 : END DO
136 52025 : DEALLOCATE (sap_int(i)%alist)
137 : END IF
138 85526 : IF (ASSOCIATED(sap_int(i)%asort)) THEN
139 51242 : DEALLOCATE (sap_int(i)%asort)
140 : END IF
141 107608 : IF (ASSOCIATED(sap_int(i)%aindex)) THEN
142 51242 : DEALLOCATE (sap_int(i)%aindex)
143 : END IF
144 : END DO
145 :
146 22082 : DEALLOCATE (sap_int)
147 :
148 22082 : 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 11020480 : 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 11020480 : NULLIFY (alist)
167 11020480 : i = locate(sap_int%asort, atom)
168 11020480 : IF (i > 0 .AND. i <= SIZE(sap_int%alist)) THEN
169 11019419 : i = sap_int%aindex(i)
170 11019419 : 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 11020480 : 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 8883074 : 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 8883074 : inn = MOD(in, 4)
201 8883074 : inn1 = inn + 1
202 63837517 : DO j = 1, jn
203 54954443 : j0 = jlist(j)
204 88579619 : DO i = 1, inn
205 33625176 : i0 = ilist(i)
206 88579619 : blk_out(i, j) = blk_in(i0, j0)
207 : END DO
208 63837517 : DO i = inn1, in, 4
209 97681786 : i0 = ilist(i)
210 97681786 : i1 = ilist(i + 1)
211 97681786 : i2 = ilist(i + 2)
212 97681786 : i3 = ilist(i + 3)
213 97681786 : blk_out(i, j) = blk_in(i0, j0)
214 97681786 : blk_out(i + 1, j) = blk_in(i1, j0)
215 97681786 : blk_out(i + 2, j) = blk_in(i2, j0)
216 104030035 : blk_out(i + 3, j) = blk_in(i3, j0)
217 : END DO
218 : END DO
219 8883074 : 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 8630672 : 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 8630672 : inn = MOD(in, 4)
242 8630672 : inn1 = inn + 1
243 61766353 : DO j = 1, jn
244 53135681 : j0 = jlist(j)
245 86532660 : DO i = 1, inn
246 33396979 : i0 = ilist(i)
247 86532660 : blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j)
248 : END DO
249 61766353 : DO i = inn1, in, 4
250 93371466 : i0 = ilist(i)
251 93371466 : i1 = ilist(i + 1)
252 93371466 : i2 = ilist(i + 2)
253 93371466 : i3 = ilist(i + 3)
254 93371466 : blk_out(i0, j0) = blk_out(i0, j0) + blk_in(i, j)
255 93371466 : blk_out(i1, j0) = blk_out(i1, j0) + blk_in(i + 1, j)
256 93371466 : blk_out(i2, j0) = blk_out(i2, j0) + blk_in(i + 2, j)
257 99298994 : blk_out(i3, j0) = blk_out(i3, j0) + blk_in(i + 3, j)
258 : END DO
259 : END DO
260 8630672 : 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 21708 : 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 21708 : !$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 21708 : 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) :: rac, raf, 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)
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 : raf(:) = pbc(particle_set(iatom)%r(:) - rf, cell)
547 : ELSE
548 : raf(:) = -rac
549 : END IF
550 : ! Keep both moment centers on the image used by the projector neighbor list.
551 : rcf(:) = raf + rac
552 :
553 : dac = NORM2(rac)
554 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
555 : clist%catom = katom
556 : clist%cell = cell_c
557 : clist%rac = rac
558 : ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
559 : clist%achint(nsgfa, nppnl, maxder))
560 : clist%acint = 0._dp
561 : clist%achint = 0._dp
562 : clist%nsgf_cnt = 0
563 : NULLIFY (clist%sgf_list)
564 :
565 : IF (PRESENT(pseudoatom)) THEN
566 : IF (pseudoatom /= katom) THEN
567 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
568 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
569 : IF (.NOT. dogth .AND. .NOT. do_upf) DEALLOCATE (radp)
570 : CYCLE
571 : END IF
572 : END IF
573 :
574 : DO iset = 1, nseta
575 : ncoa = npgfa(iset)*ncoset(la_max(iset))
576 : sgfa = first_sgfa(1, iset)
577 : IF (dogth) THEN
578 : ! GTH potential
579 : prjc = 1
580 : work = 0._dp
581 : DO l = 0, lppnl
582 : nprjc = nprj_ppnl(l)*nco(l)
583 : IF (nprjc == 0) CYCLE
584 : rprjc(1) = ppnl_radius
585 : IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
586 : lc_max = l + 2*(nprj_ppnl(l) - 1)
587 : lc_min = l
588 : zetc(1) = alpha_ppnl(l)
589 : ncoc = ncoset(lc_max)
590 :
591 : ! *** Calculate the primitive overlap integrals ***
592 : IF (my_momentmode) THEN
593 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
594 : lc_max, 1, zetc, rprjc, nder, raf, rcf, ai_work, overlap=sab)
595 : ! reduce ai_work to sab
596 : na = ncoa
597 : np = ncoc
598 : DO i = 1, maxder - 1
599 : first_col = i*ldsab
600 : sab(1:na, first_col + 1:first_col + np) = ai_work(1:na, 1:np, i)
601 : END DO
602 : ELSE
603 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
604 : lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .TRUE., ai_work, ldai)
605 : END IF
606 :
607 : ! *** Transformation step projector functions (cartesian->spherical) ***
608 : na = ncoa
609 : nb = nprjc
610 : np = ncoc
611 : DO i = 1, maxder
612 : first_col = (i - 1)*ldsab
613 : work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
614 : MATMUL(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
615 : END DO
616 : prjc = prjc + nprjc
617 : END DO
618 :
619 : na = nsgf_seta(iset)
620 : nb = nppnl
621 : np = ncoa
622 : DO i = 1, maxder
623 : first_col = (i - 1)*ldsab + 1
624 :
625 : ! *** Contraction step (basis functions) ***
626 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
627 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))
628 :
629 : ! *** Multiply with interaction matrix(h) ***
630 : clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
631 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))
632 : END DO
633 : ELSE IF (do_upf) THEN
634 : na = nsgf_seta(iset)
635 : nb = nppnl
636 : np = ncoa
637 : ALLOCATE (upf_int(np, nb, maxder))
638 : IF (my_momentmode) THEN
639 : CALL upf_projector_overlap(upf, la_max(iset), zeta(:npgfa(iset), iset), &
640 : rac, nder, upf_int, moment_origin=raf)
641 : ELSE
642 : CALL upf_projector_overlap(upf, la_max(iset), zeta(:npgfa(iset), iset), &
643 : rac, nder, upf_int)
644 : END IF
645 : DO i = 1, maxder
646 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
647 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), upf_int(:, :, i))
648 : clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
649 : MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), upf%d)
650 : END DO
651 : DEALLOCATE (upf_int)
652 : ELSE
653 : ! SGP potential
654 : ! *** Calculate the primitive overlap integrals ***
655 : IF (my_momentmode) THEN
656 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
657 : lppnl, nnl, a_nl, radp, nder, raf, rcf, ai_work, overlap=sab)
658 : ! reduce ai_work to sab
659 : na = ncoa
660 : DO i = 1, maxder - 1
661 : first_col = i*ldsab + 1
662 : sab(1:na, first_col:first_col + nprjc - 1) = ai_work(1:na, 1:nprjc, i)
663 : END DO
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 :
669 : na = nsgf_seta(iset)
670 : nb = nppnl
671 : np = ncoa
672 : DO i = 1, maxder
673 : first_col = (i - 1)*ldsab + 1
674 : ! *** Transformation step projector functions (cartesian->spherical) ***
675 : work(1:np, 1:nb) = MATMUL(sab(1:np, first_col:first_col + nprjc - 1), cprj(1:nprjc, 1:nb))
676 :
677 : ! *** Contraction step (basis functions) ***
678 : clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
679 : MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, 1:nb))
680 :
681 : ! *** Multiply with interaction matrix(h) ***
682 : ncoc = sgfa + nsgf_seta(iset) - 1
683 : DO j = 1, nppnl
684 : clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
685 : END DO
686 : END DO
687 : END IF
688 : END DO
689 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
690 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
691 : IF (.NOT. dogth .AND. .NOT. do_upf) DEALLOCATE (radp)
692 :
693 : END DO
694 :
695 : DEALLOCATE (sab, ai_work, work)
696 :
697 : !$OMP END PARALLEL
698 :
699 134 : DEALLOCATE (basis_set, gpotential, spotential)
700 :
701 134 : CALL timestop(handle)
702 :
703 402 : END SUBROUTINE build_sap_ints
704 :
705 0 : END MODULE sap_kind_types
|
