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 Calculation of the local pseudopotential contribution to the core Hamiltonian
9 : !> <a|V(local)|b> = <a|Sum e^a*rc**2|b>
10 : !> \par History
11 : !> - core_ppnl refactored from qs_core_hamiltonian [Joost VandeVondele, 2008-11-01]
12 : !> - adapted for PPL [jhu, 2009-02-23]
13 : !> - OpenMP added [Iain Bethune, Fiona Reid, 2013-11-13]
14 : !> - Bug fix: correct orbital pointer range [07.2014,JGH]
15 : !> - k-point aware [07.2015,JGH]
16 : !> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
17 : ! **************************************************************************************************
18 : MODULE core_ppl
19 :
20 : USE ai_overlap_ppl, ONLY: ecploc_integral,&
21 : ppl_integral,&
22 : ppl_integral_ri
23 : USE atomic_kind_types, ONLY: atomic_kind_type,&
24 : get_atomic_kind_set
25 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
26 : gto_basis_set_type
27 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
28 : dbcsr_get_block_p,&
29 : dbcsr_p_type
30 : USE external_potential_types, ONLY: get_potential,&
31 : gth_potential_type,&
32 : sgp_potential_type
33 : USE kinds, ONLY: dp,&
34 : int_8
35 : USE libgrpp_integrals, ONLY: libgrpp_local_forces_ref,&
36 : libgrpp_local_integrals,&
37 : libgrpp_semilocal_forces_ref,&
38 : libgrpp_semilocal_integrals
39 : USE lri_environment_types, ONLY: lri_kind_type
40 : USE orbital_pointers, ONLY: init_orbital_pointers,&
41 : ncoset
42 : USE particle_types, ONLY: particle_type
43 : USE qs_force_types, ONLY: qs_force_type
44 : USE qs_kind_types, ONLY: get_qs_kind,&
45 : get_qs_kind_set,&
46 : qs_kind_type
47 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
48 : neighbor_list_iterator_create,&
49 : neighbor_list_iterator_p_type,&
50 : neighbor_list_iterator_release,&
51 : neighbor_list_set_p_type,&
52 : nl_set_sub_iterator,&
53 : nl_sub_iterate
54 : USE virial_methods, ONLY: virial_pair_force
55 : USE virial_types, ONLY: virial_type
56 :
57 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
58 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
59 : !$ omp_init_lock, omp_set_lock, &
60 : !$ omp_unset_lock, omp_destroy_lock
61 :
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ppl'
69 :
70 : PUBLIC :: build_core_ppl, build_core_ppl_ri
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param matrix_h ...
77 : !> \param matrix_p ...
78 : !> \param force ...
79 : !> \param virial ...
80 : !> \param calculate_forces ...
81 : !> \param use_virial ...
82 : !> \param nder ...
83 : !> \param qs_kind_set ...
84 : !> \param atomic_kind_set ...
85 : !> \param particle_set ...
86 : !> \param sab_orb ...
87 : !> \param sac_ppl ...
88 : !> \param nimages ...
89 : !> \param cell_to_index ...
90 : !> \param basis_type ...
91 : !> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
92 : !> \param atcore ...
93 : ! **************************************************************************************************
94 19891 : SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
95 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
96 19891 : nimages, cell_to_index, basis_type, deltaR, atcore)
97 :
98 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
99 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
100 : TYPE(virial_type), POINTER :: virial
101 : LOGICAL, INTENT(IN) :: calculate_forces
102 : LOGICAL :: use_virial
103 : INTEGER :: nder
104 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
105 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
106 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
107 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
108 : POINTER :: sab_orb, sac_ppl
109 : INTEGER, INTENT(IN) :: nimages
110 : INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER :: cell_to_index
111 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
112 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
113 : OPTIONAL :: deltaR
114 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
115 : OPTIONAL :: atcore
116 :
117 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppl'
118 : INTEGER, PARAMETER :: nexp_max = 30
119 :
120 : INTEGER :: atom_a, handle, i, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, &
121 : katom, kkind, ldai, ldsab, maxco, maxder, maxl, maxlgto, maxlppl, maxnset, maxsgf, mepos, &
122 : n_local, natom, ncoa, ncob, nexp_lpot, nexp_ppl, nkind, nloc, nseta, nsetb, nthread, &
123 : sgfa, sgfb, slmax, slot
124 19891 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
125 : INTEGER, DIMENSION(0:10) :: npot
126 : INTEGER, DIMENSION(1:10) :: nrloc
127 : INTEGER, DIMENSION(1:15, 0:10) :: nrpot
128 : INTEGER, DIMENSION(3) :: cellind
129 19891 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, &
130 19891 : nct_lpot, npgfa, npgfb, nsgfa, nsgfb
131 19891 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
132 : INTEGER, DIMENSION(nexp_max) :: nct_ppl
133 : LOGICAL :: do_dR, doat, dokp, ecp_local, &
134 : ecp_semi_local, found, libgrpp_local, &
135 : lpotextended, only_gaussians
136 : REAL(KIND=dp) :: alpha, atk0, atk1, dab, dac, dbc, f0, &
137 : ppl_radius
138 19891 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
139 19891 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work
140 19891 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
141 : REAL(KIND=dp), ALLOCATABLE, &
142 19891 : DIMENSION(:, :, :, :, :) :: hab2
143 : REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
144 : REAL(KIND=dp), DIMENSION(1:15, 0:10) :: apot, bpot
145 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
146 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
147 : TYPE(neighbor_list_iterator_p_type), &
148 19891 : DIMENSION(:), POINTER :: ap_iterator
149 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
150 19891 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
151 : TYPE(gth_potential_type), POINTER :: gth_potential
152 39782 : REAL(KIND=dp), DIMENSION(SIZE(particle_set)) :: at_thread
153 : REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
154 19891 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h1_1block, h1_2block, &
155 19891 : h1_3block, h_block, p_block, rpgfa, &
156 19891 : rpgfb, sphi_a, sphi_b, zeta, zetb
157 19891 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
158 19891 : set_radius_a, set_radius_b
159 : REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
160 39782 : REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
161 : TYPE(sgp_potential_type), POINTER :: sgp_potential
162 :
163 : !$ INTEGER(kind=omp_lock_kind), &
164 19891 : !$ ALLOCATABLE, DIMENSION(:) :: locks
165 : !$ INTEGER :: lock_num, hash, hash1, hash2
166 : !$ INTEGER(KIND=int_8) :: iatom8
167 : !$ INTEGER, PARAMETER :: nlock = 501
168 :
169 19891 : do_dR = PRESENT(deltaR)
170 19891 : doat = PRESENT(atcore)
171 19891 : IF ((calculate_forces .OR. doat) .AND. do_dR) THEN
172 0 : CPABORT("core_ppl: incompatible options")
173 : END IF
174 :
175 : MARK_USED(int_8)
176 :
177 : ! Use internal integral routine for local ECP terms or use libgrrp
178 19891 : libgrpp_local = .FALSE.
179 :
180 19891 : IF (calculate_forces) THEN
181 7919 : CALL timeset(routineN//"_forces", handle)
182 : ELSE
183 11972 : CALL timeset(routineN, handle)
184 : END IF
185 :
186 19891 : nkind = SIZE(atomic_kind_set)
187 19891 : natom = SIZE(particle_set)
188 :
189 19891 : dokp = (nimages > 1)
190 :
191 19891 : IF (dokp) THEN
192 658 : IF (PRESENT(cell_to_index)) THEN
193 658 : CPASSERT(ASSOCIATED(cell_to_index))
194 : ELSE
195 0 : CPABORT("Missing cell_to_index for k-point calculation")
196 : END IF
197 : END IF
198 :
199 19891 : IF (calculate_forces .OR. doat) THEN
200 7981 : IF (SIZE(matrix_p, 1) == 2) THEN
201 2550 : DO img = 1, nimages
202 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
203 1576 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
204 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
205 2550 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
206 : END DO
207 : END IF
208 : END IF
209 318595 : force_thread = 0.0_dp
210 94567 : at_thread = 0.0_dp
211 :
212 19891 : maxder = ncoset(nder)
213 :
214 : CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxlgto=maxlgto, &
215 : maxsgf=maxsgf, maxnset=maxnset, maxlppl=maxlppl, &
216 19891 : basis_type=basis_type)
217 :
218 19891 : maxl = MAX(maxlgto, maxlppl)
219 19891 : CALL init_orbital_pointers(2*maxl + 2*nder + 1)
220 :
221 19891 : ldsab = MAX(maxco, ncoset(maxlppl), maxsgf, maxlppl)
222 19891 : ldai = ncoset(maxl + nder + 1)
223 :
224 94732 : ALLOCATE (basis_set_list(nkind))
225 54950 : DO ikind = 1, nkind
226 35059 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
227 54950 : IF (ASSOCIATED(basis_set_a)) THEN
228 35059 : basis_set_list(ikind)%gto_basis_set => basis_set_a
229 : ELSE
230 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
231 : END IF
232 : END DO
233 :
234 19891 : pv_thread = 0.0_dp
235 :
236 : nthread = 1
237 19891 : !$ nthread = omp_get_max_threads()
238 :
239 : ! iterator for basis/potential list
240 19891 : CALL neighbor_list_iterator_create(ap_iterator, sac_ppl, search=.TRUE., nthread=nthread)
241 :
242 : !$OMP PARALLEL &
243 : !$OMP DEFAULT (NONE) &
244 : !$OMP SHARED (ap_iterator, basis_set_list, calculate_forces, use_virial, &
245 : !$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
246 : !$OMP sab_orb, sac_ppl, nthread, ncoset, nkind, cell_to_index, &
247 : !$OMP ldsab, maxnset, maxder, do_dR, deltaR, doat, libgrpp_local, &
248 : !$OMP maxlgto, nder, maxco, dokp, locks, natom) &
249 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
250 : !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
251 : !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
252 : !$OMP zetb, dab, irow, icol, h_block, found, iset, ncoa, lock_num, &
253 : !$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, hab2, hab2_w, &
254 : !$OMP atk0, atk1, h1_1block, h1_2block, h1_3block, kkind, nseta, &
255 : !$OMP gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended, &
256 : !$OMP ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl, &
257 : !$OMP nct_lpot, cval_ppl, cval_lpot, rac, dac, rbc, dbc, &
258 : !$OMP set_radius_a, rpgfa, force_a, force_b, ppl_fwork, mepos, &
259 : !$OMP slot, f0, katom, ppl_work, cellind, img, ecp_local, ecp_semi_local, &
260 : !$OMP nloc, nrloc, aloc, bloc, n_local, a_local, c_local, &
261 : !$OMP slmax, npot, nrpot, apot, bpot, only_gaussians, &
262 : !$OMP ldai, hash, hash1, hash2, iatom8) &
263 19891 : !$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )
264 :
265 : !$OMP SINGLE
266 : !$ ALLOCATE (locks(nlock))
267 : !$OMP END SINGLE
268 :
269 : !$OMP DO
270 : !$ DO lock_num = 1, nlock
271 : !$ call omp_init_lock(locks(lock_num))
272 : !$ END DO
273 : !$OMP END DO
274 :
275 : mepos = 0
276 : !$ mepos = omp_get_thread_num()
277 :
278 : ALLOCATE (hab(ldsab, ldsab, maxnset, maxnset), work(ldsab, ldsab*maxder))
279 : ldai = ncoset(2*maxlgto + 2*nder)
280 : ALLOCATE (ppl_work(ldai, ldai, MAX(maxder, 2*maxlgto + 2*nder + 1)))
281 : IF (calculate_forces .OR. doat) THEN
282 : ALLOCATE (pab(maxco, maxco, maxnset, maxnset))
283 : ldai = ncoset(maxlgto)
284 : ALLOCATE (ppl_fwork(ldai, ldai, maxder))
285 : END IF
286 :
287 : !$OMP DO SCHEDULE(GUIDED)
288 : DO slot = 1, sab_orb(1)%nl_size
289 : !SL
290 : IF (do_dR) THEN
291 : ALLOCATE (hab2(ldsab, ldsab, 4, maxnset, maxnset))
292 : ALLOCATE (hab2_w(ldsab, ldsab, 6))
293 : ALLOCATE (ppl_fwork(ldai, ldai, maxder))
294 : END IF
295 :
296 : ikind = sab_orb(1)%nlist_task(slot)%ikind
297 : jkind = sab_orb(1)%nlist_task(slot)%jkind
298 : iatom = sab_orb(1)%nlist_task(slot)%iatom
299 : jatom = sab_orb(1)%nlist_task(slot)%jatom
300 : cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
301 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
302 :
303 : basis_set_a => basis_set_list(ikind)%gto_basis_set
304 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
305 : basis_set_b => basis_set_list(jkind)%gto_basis_set
306 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
307 :
308 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
309 : !$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
310 :
311 : ! basis ikind
312 : first_sgfa => basis_set_a%first_sgf
313 : la_max => basis_set_a%lmax
314 : la_min => basis_set_a%lmin
315 : npgfa => basis_set_a%npgf
316 : nseta = basis_set_a%nset
317 : nsgfa => basis_set_a%nsgf_set
318 : rpgfa => basis_set_a%pgf_radius
319 : set_radius_a => basis_set_a%set_radius
320 : sphi_a => basis_set_a%sphi
321 : zeta => basis_set_a%zet
322 : ! basis jkind
323 : first_sgfb => basis_set_b%first_sgf
324 : lb_max => basis_set_b%lmax
325 : lb_min => basis_set_b%lmin
326 : npgfb => basis_set_b%npgf
327 : nsetb = basis_set_b%nset
328 : nsgfb => basis_set_b%nsgf_set
329 : rpgfb => basis_set_b%pgf_radius
330 : set_radius_b => basis_set_b%set_radius
331 : sphi_b => basis_set_b%sphi
332 : zetb => basis_set_b%zet
333 :
334 : dab = SQRT(SUM(rab*rab))
335 :
336 : IF (dokp) THEN
337 : img = cell_to_index(cellind(1), cellind(2), cellind(3))
338 : ELSE
339 : img = 1
340 : END IF
341 :
342 : ! *** Use the symmetry of the first derivatives ***
343 : IF (iatom == jatom) THEN
344 : f0 = 1.0_dp
345 : ELSE
346 : f0 = 2.0_dp
347 : END IF
348 :
349 : ! *** Create matrix blocks for a new matrix block column ***
350 : IF (iatom <= jatom) THEN
351 : irow = iatom
352 : icol = jatom
353 : ELSE
354 : irow = jatom
355 : icol = iatom
356 : END IF
357 : NULLIFY (h_block)
358 :
359 : IF (do_dR) THEN
360 : NULLIFY (h1_1block, h1_2block, h1_3block)
361 :
362 : CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, &
363 : row=irow, col=icol, BLOCK=h1_1block, found=found)
364 : CALL dbcsr_get_block_p(matrix=matrix_h(2, img)%matrix, &
365 : row=irow, col=icol, BLOCK=h1_2block, found=found)
366 : CALL dbcsr_get_block_p(matrix=matrix_h(3, img)%matrix, &
367 : row=irow, col=icol, BLOCK=h1_3block, found=found)
368 : END IF
369 :
370 : CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
371 : CPASSERT(found)
372 : IF (calculate_forces .OR. doat) THEN
373 : NULLIFY (p_block)
374 : CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
375 : IF (ASSOCIATED(p_block)) THEN
376 : DO iset = 1, nseta
377 : ncoa = npgfa(iset)*ncoset(la_max(iset))
378 : sgfa = first_sgfa(1, iset)
379 : DO jset = 1, nsetb
380 : ncob = npgfb(jset)*ncoset(lb_max(jset))
381 : sgfb = first_sgfb(1, jset)
382 :
383 : ! *** Decontract density matrix block ***
384 : IF (iatom <= jatom) THEN
385 : work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
386 : p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
387 : ELSE
388 : work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
389 : TRANSPOSE(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1)))
390 : END IF
391 :
392 : pab(1:ncoa, 1:ncob, iset, jset) = MATMUL(work(1:ncoa, 1:nsgfb(jset)), &
393 : TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
394 : END DO
395 : END DO
396 : END IF
397 : END IF
398 :
399 : hab = 0._dp
400 : IF (do_dr) hab2 = 0._dp
401 :
402 : ! loop over all kinds for pseudopotential atoms
403 : DO kkind = 1, nkind
404 :
405 : CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
406 : sgp_potential=sgp_potential)
407 : ecp_semi_local = .FALSE.
408 : only_gaussians = .TRUE.
409 : IF (ASSOCIATED(gth_potential)) THEN
410 : CALL get_potential(potential=gth_potential, &
411 : alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
412 : lpot_present=lpotextended, ppl_radius=ppl_radius)
413 : nexp_ppl = 1
414 : alpha_ppl(1) = alpha
415 : nct_ppl(1) = SIZE(cexp_ppl)
416 : cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
417 : IF (lpotextended) THEN
418 : CALL get_potential(potential=gth_potential, &
419 : nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, &
420 : cval_lpot=cval_lpot)
421 : CPASSERT(nexp_lpot < nexp_max)
422 : nexp_ppl = nexp_lpot + 1
423 : alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
424 : nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
425 : DO i = 1, nexp_lpot
426 : cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
427 : END DO
428 : END IF
429 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
430 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
431 : ppl_radius=ppl_radius)
432 : IF (ecp_local) THEN
433 : CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
434 : nexp_ppl = nloc
435 : CPASSERT(nexp_ppl <= nexp_max)
436 : nct_ppl(1:nloc) = nrloc(1:nloc)
437 : alpha_ppl(1:nloc) = bloc(1:nloc)
438 : cval_ppl(1, 1:nloc) = aloc(1:nloc)
439 : only_gaussians = .FALSE.
440 : ELSE
441 : CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
442 : nexp_ppl = n_local
443 : CPASSERT(nexp_ppl <= nexp_max)
444 : nct_ppl(1:n_local) = 1
445 : alpha_ppl(1:n_local) = a_local(1:n_local)
446 : cval_ppl(1, 1:n_local) = c_local(1:n_local)
447 : END IF
448 : IF (ecp_semi_local) THEN
449 : CALL get_potential(potential=sgp_potential, sl_lmax=slmax, &
450 : npot=npot, nrpot=nrpot, apot=apot, bpot=bpot)
451 : ELSE IF (ecp_local) THEN
452 : IF (SUM(ABS(aloc(1:nloc))) < 1.0e-12_dp) CYCLE
453 : END IF
454 : ELSE
455 : CYCLE
456 : END IF
457 :
458 : CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
459 :
460 : DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
461 :
462 : CALL get_iterator_info(ap_iterator, mepos=mepos, jatom=katom, r=rac)
463 :
464 : dac = SQRT(SUM(rac*rac))
465 : rbc(:) = rac(:) - rab(:)
466 : dbc = SQRT(SUM(rbc*rbc))
467 : IF ((MAXVAL(set_radius_a(:)) + ppl_radius < dac) .OR. &
468 : (MAXVAL(set_radius_b(:)) + ppl_radius < dbc)) THEN
469 : CYCLE
470 : END IF
471 :
472 : DO iset = 1, nseta
473 : IF (set_radius_a(iset) + ppl_radius < dac) CYCLE
474 : ncoa = npgfa(iset)*ncoset(la_max(iset))
475 : sgfa = first_sgfa(1, iset)
476 : DO jset = 1, nsetb
477 : IF (set_radius_b(jset) + ppl_radius < dbc) CYCLE
478 : ncob = npgfb(jset)*ncoset(lb_max(jset))
479 : sgfb = first_sgfb(1, jset)
480 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
481 : ! *** Calculate the GTH pseudo potential forces ***
482 : IF (doat) THEN
483 : atk0 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
484 : pab(1:ncoa, 1:ncob, iset, jset))
485 : END IF
486 : IF (calculate_forces) THEN
487 :
488 : force_a(:) = 0.0_dp
489 : force_b(:) = 0.0_dp
490 :
491 : IF (only_gaussians) THEN
492 : CALL ppl_integral( &
493 : la_max(iset), la_min(iset), npgfa(iset), &
494 : rpgfa(:, iset), zeta(:, iset), &
495 : lb_max(jset), lb_min(jset), npgfb(jset), &
496 : rpgfb(:, jset), zetb(:, jset), &
497 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
498 : rab, dab, rac, dac, rbc, dbc, &
499 : hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
500 : force_a, force_b, ppl_fwork)
501 : ELSE IF (libgrpp_local) THEN
502 : !$OMP CRITICAL(type1)
503 : CALL libgrpp_local_forces_ref(la_max(iset), la_min(iset), npgfa(iset), &
504 : rpgfa(:, iset), zeta(:, iset), &
505 : lb_max(jset), lb_min(jset), npgfb(jset), &
506 : rpgfb(:, jset), zetb(:, jset), &
507 : nexp_ppl, alpha_ppl, cval_ppl(1, :), nct_ppl, &
508 : ppl_radius, rab, dab, rac, dac, dbc, &
509 : hab(:, :, iset, jset), pab(:, :, iset, jset), &
510 : force_a, force_b)
511 : !$OMP END CRITICAL(type1)
512 : ELSE
513 : CALL ecploc_integral( &
514 : la_max(iset), la_min(iset), npgfa(iset), &
515 : rpgfa(:, iset), zeta(:, iset), &
516 : lb_max(jset), lb_min(jset), npgfb(jset), &
517 : rpgfb(:, jset), zetb(:, jset), &
518 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
519 : rab, dab, rac, dac, rbc, dbc, &
520 : hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
521 : force_a, force_b, ppl_fwork)
522 : END IF
523 :
524 : IF (ecp_semi_local) THEN
525 :
526 : !$OMP CRITICAL(type2)
527 : CALL libgrpp_semilocal_forces_ref(la_max(iset), la_min(iset), npgfa(iset), &
528 : rpgfa(:, iset), zeta(:, iset), &
529 : lb_max(jset), lb_min(jset), npgfb(jset), &
530 : rpgfb(:, jset), zetb(:, jset), &
531 : slmax, npot, bpot, apot, nrpot, &
532 : ppl_radius, rab, dab, rac, dac, dbc, &
533 : hab(:, :, iset, jset), pab(:, :, iset, jset), &
534 : force_a, force_b)
535 : !$OMP END CRITICAL(type2)
536 : END IF
537 : ! *** The derivatives w.r.t. atomic center c are ***
538 : ! *** calculated using the translational invariance ***
539 : ! *** of the first derivatives ***
540 :
541 : force_thread(1, iatom) = force_thread(1, iatom) + f0*force_a(1)
542 : force_thread(2, iatom) = force_thread(2, iatom) + f0*force_a(2)
543 : force_thread(3, iatom) = force_thread(3, iatom) + f0*force_a(3)
544 : force_thread(1, katom) = force_thread(1, katom) - f0*force_a(1)
545 : force_thread(2, katom) = force_thread(2, katom) - f0*force_a(2)
546 : force_thread(3, katom) = force_thread(3, katom) - f0*force_a(3)
547 :
548 : force_thread(1, jatom) = force_thread(1, jatom) + f0*force_b(1)
549 : force_thread(2, jatom) = force_thread(2, jatom) + f0*force_b(2)
550 : force_thread(3, jatom) = force_thread(3, jatom) + f0*force_b(3)
551 : force_thread(1, katom) = force_thread(1, katom) - f0*force_b(1)
552 : force_thread(2, katom) = force_thread(2, katom) - f0*force_b(2)
553 : force_thread(3, katom) = force_thread(3, katom) - f0*force_b(3)
554 :
555 : IF (use_virial) THEN
556 : CALL virial_pair_force(pv_thread, f0, force_a, rac)
557 : CALL virial_pair_force(pv_thread, f0, force_b, rbc)
558 : END IF
559 : ELSE IF (do_dR) THEN
560 : hab2_w = 0._dp
561 : CALL ppl_integral( &
562 : la_max(iset), la_min(iset), npgfa(iset), &
563 : rpgfa(:, iset), zeta(:, iset), &
564 : lb_max(jset), lb_min(jset), npgfb(jset), &
565 : rpgfb(:, jset), zetb(:, jset), &
566 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
567 : rab, dab, rac, dac, rbc, dbc, &
568 : vab=hab(:, :, iset, jset), s=ppl_work, &
569 : hab2=hab2(:, :, :, iset, jset), hab2_work=hab2_w, fs=ppl_fwork, &
570 : deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
571 : IF (ecp_semi_local) THEN
572 : ! semi local ECP part
573 : CPABORT("Option not implemented")
574 : END IF
575 : ELSE
576 : IF (only_gaussians) THEN
577 : !If the local part of the pseudo-potential only has Gaussian functions
578 : !we can use CP2K native code, that can run without libgrpp installation
579 : CALL ppl_integral( &
580 : la_max(iset), la_min(iset), npgfa(iset), &
581 : rpgfa(:, iset), zeta(:, iset), &
582 : lb_max(jset), lb_min(jset), npgfb(jset), &
583 : rpgfb(:, jset), zetb(:, jset), &
584 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
585 : rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
586 :
587 : ELSE IF (libgrpp_local) THEN
588 : !If the local part of the potential is more complex, we need libgrpp
589 : !$OMP CRITICAL(type1)
590 : CALL libgrpp_local_integrals(la_max(iset), la_min(iset), npgfa(iset), &
591 : rpgfa(:, iset), zeta(:, iset), &
592 : lb_max(jset), lb_min(jset), npgfb(jset), &
593 : rpgfb(:, jset), zetb(:, jset), &
594 : nexp_ppl, alpha_ppl, cval_ppl(1, :), nct_ppl, &
595 : ppl_radius, rab, dab, rac, dac, dbc, &
596 : hab(:, :, iset, jset))
597 : !$OMP END CRITICAL(type1)
598 : ELSE
599 : CALL ecploc_integral( &
600 : la_max(iset), la_min(iset), npgfa(iset), &
601 : rpgfa(:, iset), zeta(:, iset), &
602 : lb_max(jset), lb_min(jset), npgfb(jset), &
603 : rpgfb(:, jset), zetb(:, jset), &
604 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
605 : rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
606 : END IF
607 :
608 : IF (ecp_semi_local) THEN
609 : ! semi local ECP part
610 : !$OMP CRITICAL(type2)
611 : CALL libgrpp_semilocal_integrals(la_max(iset), la_min(iset), npgfa(iset), &
612 : rpgfa(:, iset), zeta(:, iset), &
613 : lb_max(jset), lb_min(jset), npgfb(jset), &
614 : rpgfb(:, jset), zetb(:, jset), &
615 : slmax, npot, bpot, apot, nrpot, &
616 : ppl_radius, rab, dab, rac, dac, dbc, &
617 : hab(:, :, iset, jset))
618 : !$OMP END CRITICAL(type2)
619 : END IF
620 : END IF
621 : ! calculate atomic contributions
622 : IF (doat) THEN
623 : atk1 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
624 : pab(1:ncoa, 1:ncob, iset, jset))
625 : at_thread(katom) = at_thread(katom) + (atk1 - atk0)
626 : END IF
627 : END DO
628 : END DO
629 : END DO
630 : END DO
631 :
632 : ! *** Contract PPL integrals
633 : IF (.NOT. do_dR) THEN
634 : DO iset = 1, nseta
635 : ncoa = npgfa(iset)*ncoset(la_max(iset))
636 : sgfa = first_sgfa(1, iset)
637 : DO jset = 1, nsetb
638 : ncob = npgfb(jset)*ncoset(lb_max(jset))
639 : sgfb = first_sgfb(1, jset)
640 :
641 : !$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
642 : !$ hash = MOD(hash1 + hash2, nlock) + 1
643 :
644 : work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab(1:ncoa, 1:ncob, iset, jset), &
645 : sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
646 : !$ CALL omp_set_lock(locks(hash))
647 : IF (iatom <= jatom) THEN
648 : h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
649 : h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
650 : MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
651 : ELSE
652 : h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
653 : h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
654 : MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
655 : END IF
656 : !$ CALL omp_unset_lock(locks(hash))
657 :
658 : END DO
659 : END DO
660 : ELSE ! do_dr == .true.
661 : DO iset = 1, nseta
662 : ncoa = npgfa(iset)*ncoset(la_max(iset))
663 : sgfa = first_sgfa(1, iset)
664 : DO jset = 1, nsetb
665 : ncob = npgfb(jset)*ncoset(lb_max(jset))
666 : sgfb = first_sgfb(1, jset)
667 : work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 1, iset, jset), &
668 : sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
669 :
670 : !$OMP CRITICAL(h1_1block_critical)
671 : IF (iatom <= jatom) THEN
672 : h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
673 : h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
674 : MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
675 :
676 : ELSE
677 : h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
678 : h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
679 : MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
680 : END IF
681 : !$OMP END CRITICAL(h1_1block_critical)
682 : work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 2, iset, jset), &
683 : sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
684 :
685 : !$OMP CRITICAL(h1_2block_critical)
686 : IF (iatom <= jatom) THEN
687 : h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
688 : h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
689 : MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
690 :
691 : ELSE
692 : h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
693 : h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
694 : MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
695 : END IF
696 : !$OMP END CRITICAL(h1_2block_critical)
697 : work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 3, iset, jset), &
698 : sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
699 : !$OMP CRITICAL(h1_3block_critical)
700 : IF (iatom <= jatom) THEN
701 : h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
702 : h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
703 : MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
704 :
705 : ELSE
706 : h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
707 : h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
708 : MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
709 : END IF
710 : !$OMP END CRITICAL(h1_3block_critical)
711 : END DO
712 : END DO
713 : END IF
714 : IF (do_dR) DEALLOCATE (hab2, ppl_fwork, hab2_w)
715 : END DO ! slot
716 :
717 : DEALLOCATE (hab, work, ppl_work)
718 : IF (calculate_forces .OR. doat) THEN
719 : DEALLOCATE (pab, ppl_fwork)
720 : END IF
721 :
722 : !$OMP DO
723 : !$ DO lock_num = 1, nlock
724 : !$ call omp_destroy_lock(locks(lock_num))
725 : !$ END DO
726 : !$OMP END DO
727 :
728 : !$OMP SINGLE
729 : !$ DEALLOCATE (locks)
730 : !$OMP END SINGLE NOWAIT
731 :
732 : !$OMP END PARALLEL
733 :
734 19891 : CALL neighbor_list_iterator_release(ap_iterator)
735 :
736 19891 : DEALLOCATE (basis_set_list)
737 :
738 19891 : IF (calculate_forces .OR. doat) THEN
739 : ! *** If LSD, then recover alpha density and beta density ***
740 : ! *** from the total density (1) and the spin density (2) ***
741 7981 : IF (SIZE(matrix_p, 1) == 2) THEN
742 2550 : DO img = 1, nimages
743 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
744 1576 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
745 : CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
746 2550 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
747 : END DO
748 : END IF
749 : END IF
750 :
751 19891 : IF (calculate_forces) THEN
752 7919 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
753 : !$OMP DO
754 : DO iatom = 1, natom
755 28231 : atom_a = atom_of_kind(iatom)
756 28231 : ikind = kind_of(iatom)
757 112924 : force(ikind)%gth_ppl(:, atom_a) = force(ikind)%gth_ppl(:, atom_a) + force_thread(:, iatom)
758 : END DO
759 : !$OMP END DO
760 7919 : DEALLOCATE (atom_of_kind, kind_of)
761 : END IF
762 19891 : IF (doat) THEN
763 280 : atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
764 : END IF
765 :
766 19891 : IF (calculate_forces .AND. use_virial) THEN
767 11986 : virial%pv_ppl = virial%pv_ppl + pv_thread
768 11986 : virial%pv_virial = virial%pv_virial + pv_thread
769 : END IF
770 :
771 19891 : CALL timestop(handle)
772 :
773 59673 : END SUBROUTINE build_core_ppl
774 :
775 : ! **************************************************************************************************
776 : !> \brief ...
777 : !> \param lri_ppl_coef ...
778 : !> \param force ...
779 : !> \param virial ...
780 : !> \param calculate_forces ...
781 : !> \param use_virial ...
782 : !> \param qs_kind_set ...
783 : !> \param atomic_kind_set ...
784 : !> \param particle_set ...
785 : !> \param sac_ppl ...
786 : !> \param basis_type ...
787 : ! **************************************************************************************************
788 4 : SUBROUTINE build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_virial, &
789 : qs_kind_set, atomic_kind_set, particle_set, sac_ppl, &
790 : basis_type)
791 :
792 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_ppl_coef
793 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
794 : TYPE(virial_type), POINTER :: virial
795 : LOGICAL, INTENT(IN) :: calculate_forces
796 : LOGICAL :: use_virial
797 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
798 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
799 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
800 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
801 : POINTER :: sac_ppl
802 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
803 :
804 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppl_ri'
805 : INTEGER, PARAMETER :: nexp_max = 30
806 :
807 : INTEGER :: atom_a, handle, i, iatom, ikind, iset, katom, kkind, maxco, maxsgf, n_local, &
808 : natom, ncoa, nexp_lpot, nexp_ppl, nfun, nkind, nloc, nseta, sgfa, sgfb, slot
809 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
810 : INTEGER, DIMENSION(1:10) :: nrloc
811 4 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, nct_lpot, npgfa, nsgfa
812 4 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
813 : INTEGER, DIMENSION(nexp_max) :: nct_ppl
814 : LOGICAL :: ecp_local, ecp_semi_local, lpotextended
815 : REAL(KIND=dp) :: alpha, dac, ppl_radius
816 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: va, work
817 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dva, dvas
818 : REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
819 : REAL(KIND=dp), DIMENSION(3) :: force_a, rac
820 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
821 : TYPE(gto_basis_set_type), POINTER :: basis_set
822 4 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
823 : TYPE(gth_potential_type), POINTER :: gth_potential
824 : REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
825 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: bcon, cval_lpot, rpgfa, sphi_a, zeta
826 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
827 4 : set_radius_a
828 : REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
829 8 : REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
830 : TYPE(sgp_potential_type), POINTER :: sgp_potential
831 :
832 : !$ INTEGER(kind=omp_lock_kind), &
833 4 : !$ ALLOCATABLE, DIMENSION(:) :: locks
834 : !$ INTEGER :: lock_num, hash
835 : !$ INTEGER, PARAMETER :: nlock = 501
836 :
837 8 : IF (calculate_forces) THEN
838 2 : CALL timeset(routineN//"_forces", handle)
839 : ELSE
840 2 : CALL timeset(routineN, handle)
841 : END IF
842 :
843 4 : nkind = SIZE(atomic_kind_set)
844 4 : natom = SIZE(particle_set)
845 :
846 52 : force_thread = 0.0_dp
847 4 : pv_thread = 0.0_dp
848 4 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
849 :
850 20 : ALLOCATE (basis_set_list(nkind))
851 12 : DO ikind = 1, nkind
852 8 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
853 12 : IF (ASSOCIATED(basis_set)) THEN
854 8 : basis_set_list(ikind)%gto_basis_set => basis_set
855 : ELSE
856 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
857 : END IF
858 : END DO
859 :
860 4 : CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxsgf=maxsgf, basis_type=basis_type)
861 :
862 : !$OMP PARALLEL &
863 : !$OMP DEFAULT (NONE) &
864 : !$OMP SHARED (maxco,maxsgf,basis_set_list,calculate_forces,lri_ppl_coef,qs_kind_set,&
865 : !$OMP locks,natom,use_virial,virial,ncoset,atom_of_kind,sac_ppl) &
866 : !$OMP PRIVATE (ikind,kkind,iatom,katom,atom_a,rac,va,dva,dvas,basis_set,slot,&
867 : !$OMP first_sgfa,la_max,la_min,npgfa,nseta,nsgfa,rpgfa,set_radius_a,lock_num,&
868 : !$OMP sphi_a,zeta,gth_potential,sgp_potential,alpha,cexp_ppl,lpotextended,ppl_radius,&
869 : !$OMP nexp_ppl,alpha_ppl,nct_ppl,cval_ppl,nloc,n_local,nrloc,a_local,aloc,bloc,c_local,nfun,work,&
870 : !$OMP hash,dac,force_a,iset,sgfa,sgfb,ncoa,bcon,cval_lpot,nct_lpot,alpha_lpot,nexp_lpot,&
871 : !$OMP ecp_local,ecp_semi_local) &
872 4 : !$OMP REDUCTION (+ : pv_thread, force_thread )
873 :
874 : !$OMP SINGLE
875 : !$ ALLOCATE (locks(nlock))
876 : !$OMP END SINGLE
877 :
878 : !$OMP DO
879 : !$ DO lock_num = 1, nlock
880 : !$ call omp_init_lock(locks(lock_num))
881 : !$ END DO
882 : !$OMP END DO
883 :
884 : ALLOCATE (va(maxco), work(maxsgf))
885 : IF (calculate_forces) THEN
886 : ALLOCATE (dva(maxco, 3), dvas(maxco, 3))
887 : END IF
888 :
889 : !$OMP DO SCHEDULE(GUIDED)
890 : DO slot = 1, sac_ppl(1)%nl_size
891 :
892 : ikind = sac_ppl(1)%nlist_task(slot)%ikind
893 : kkind = sac_ppl(1)%nlist_task(slot)%jkind
894 : iatom = sac_ppl(1)%nlist_task(slot)%iatom
895 : katom = sac_ppl(1)%nlist_task(slot)%jatom
896 : rac(1:3) = sac_ppl(1)%nlist_task(slot)%r(1:3)
897 : atom_a = atom_of_kind(iatom)
898 :
899 : basis_set => basis_set_list(ikind)%gto_basis_set
900 : IF (.NOT. ASSOCIATED(basis_set)) CYCLE
901 :
902 : ! basis ikind
903 : first_sgfa => basis_set%first_sgf
904 : la_max => basis_set%lmax
905 : la_min => basis_set%lmin
906 : npgfa => basis_set%npgf
907 : nseta = basis_set%nset
908 : nsgfa => basis_set%nsgf_set
909 : nfun = basis_set%nsgf
910 : rpgfa => basis_set%pgf_radius
911 : set_radius_a => basis_set%set_radius
912 : sphi_a => basis_set%sphi
913 : zeta => basis_set%zet
914 :
915 : CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
916 : sgp_potential=sgp_potential)
917 : ecp_semi_local = .FALSE.
918 : IF (ASSOCIATED(gth_potential)) THEN
919 : CALL get_potential(potential=gth_potential, &
920 : alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
921 : lpot_present=lpotextended, ppl_radius=ppl_radius)
922 : nexp_ppl = 1
923 : alpha_ppl(1) = alpha
924 : nct_ppl(1) = SIZE(cexp_ppl)
925 : cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
926 : IF (lpotextended) THEN
927 : CALL get_potential(potential=gth_potential, &
928 : nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, cval_lpot=cval_lpot)
929 : CPASSERT(nexp_lpot < nexp_max)
930 : nexp_ppl = nexp_lpot + 1
931 : alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
932 : nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
933 : DO i = 1, nexp_lpot
934 : cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
935 : END DO
936 : END IF
937 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
938 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
939 : ppl_radius=ppl_radius)
940 : CPASSERT(.NOT. ecp_semi_local)
941 : IF (ecp_local) THEN
942 : CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
943 : IF (SUM(ABS(aloc(1:nloc))) < 1.0e-12_dp) CYCLE
944 : nexp_ppl = nloc
945 : CPASSERT(nexp_ppl <= nexp_max)
946 : nct_ppl(1:nloc) = nrloc(1:nloc)
947 : alpha_ppl(1:nloc) = bloc(1:nloc)
948 : cval_ppl(1, 1:nloc) = aloc(1:nloc)
949 : ELSE
950 : CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
951 : nexp_ppl = n_local
952 : CPASSERT(nexp_ppl <= nexp_max)
953 : nct_ppl(1:n_local) = 1
954 : alpha_ppl(1:n_local) = a_local(1:n_local)
955 : cval_ppl(1, 1:n_local) = c_local(1:n_local)
956 : END IF
957 : ELSE
958 : CYCLE
959 : END IF
960 :
961 : dac = SQRT(SUM(rac*rac))
962 : IF ((MAXVAL(set_radius_a(:)) + ppl_radius < dac)) CYCLE
963 : IF (calculate_forces) force_a = 0.0_dp
964 : work(1:nfun) = 0.0_dp
965 :
966 : DO iset = 1, nseta
967 : IF (set_radius_a(iset) + ppl_radius < dac) CYCLE
968 : ! integrals
969 : IF (calculate_forces) THEN
970 : va = 0.0_dp
971 : dva = 0.0_dp
972 : CALL ppl_integral_ri( &
973 : la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
974 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
975 : -rac, dac, va, dva)
976 : ELSE
977 : va = 0.0_dp
978 : CALL ppl_integral_ri( &
979 : la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
980 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
981 : -rac, dac, va)
982 : END IF
983 : ! contraction
984 : sgfa = first_sgfa(1, iset)
985 : sgfb = sgfa + nsgfa(iset) - 1
986 : ncoa = npgfa(iset)*ncoset(la_max(iset))
987 : bcon => sphi_a(1:ncoa, sgfa:sgfb)
988 : work(sgfa:sgfb) = MATMUL(TRANSPOSE(bcon), va(1:ncoa))
989 : IF (calculate_forces) THEN
990 : dvas(1:nsgfa(iset), 1:3) = MATMUL(TRANSPOSE(bcon), dva(1:ncoa, 1:3))
991 : force_a(1) = force_a(1) + SUM(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 1))
992 : force_a(2) = force_a(2) + SUM(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 2))
993 : force_a(3) = force_a(3) + SUM(lri_ppl_coef(ikind)%acoef(atom_a, sgfa:sgfb)*dvas(1:nsgfa(iset), 3))
994 : END IF
995 : END DO
996 : !$ hash = MOD(iatom, nlock) + 1
997 : !$ CALL omp_set_lock(locks(hash))
998 : lri_ppl_coef(ikind)%v_int(atom_a, 1:nfun) = lri_ppl_coef(ikind)%v_int(atom_a, 1:nfun) + work(1:nfun)
999 : !$ CALL omp_unset_lock(locks(hash))
1000 : IF (calculate_forces) THEN
1001 : force_thread(1, iatom) = force_thread(1, iatom) + force_a(1)
1002 : force_thread(2, iatom) = force_thread(2, iatom) + force_a(2)
1003 : force_thread(3, iatom) = force_thread(3, iatom) + force_a(3)
1004 : force_thread(1, katom) = force_thread(1, katom) - force_a(1)
1005 : force_thread(2, katom) = force_thread(2, katom) - force_a(2)
1006 : force_thread(3, katom) = force_thread(3, katom) - force_a(3)
1007 : IF (use_virial) THEN
1008 : CALL virial_pair_force(pv_thread, 1.0_dp, force_a, rac)
1009 : END IF
1010 : END IF
1011 : END DO
1012 :
1013 : DEALLOCATE (va, work)
1014 : IF (calculate_forces) THEN
1015 : DEALLOCATE (dva, dvas)
1016 : END IF
1017 :
1018 : !$OMP END PARALLEL
1019 :
1020 4 : IF (calculate_forces) THEN
1021 8 : DO iatom = 1, natom
1022 6 : atom_a = atom_of_kind(iatom)
1023 6 : ikind = kind_of(iatom)
1024 6 : force(ikind)%gth_ppl(1, atom_a) = force(ikind)%gth_ppl(1, atom_a) + force_thread(1, iatom)
1025 6 : force(ikind)%gth_ppl(2, atom_a) = force(ikind)%gth_ppl(2, atom_a) + force_thread(2, iatom)
1026 8 : force(ikind)%gth_ppl(3, atom_a) = force(ikind)%gth_ppl(3, atom_a) + force_thread(3, iatom)
1027 : END DO
1028 : END IF
1029 4 : DEALLOCATE (atom_of_kind, kind_of)
1030 :
1031 4 : IF (calculate_forces .AND. use_virial) THEN
1032 0 : virial%pv_ppl = virial%pv_ppl + pv_thread
1033 0 : virial%pv_virial = virial%pv_virial + pv_thread
1034 : END IF
1035 :
1036 4 : DEALLOCATE (basis_set_list)
1037 :
1038 4 : CALL timestop(handle)
1039 :
1040 12 : END SUBROUTINE build_core_ppl_ri
1041 :
1042 : ! **************************************************************************************************
1043 :
1044 : END MODULE core_ppl
|
