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