(git:c51e757)
Loading...
Searching...
No Matches
core_ppl.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief 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! **************************************************************************************************
18MODULE 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
72CONTAINS
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 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 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 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, &
130 nct_lpot, npgfa, npgfb, nsgfa, nsgfb
131 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
139 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work
140 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
141 REAL(kind=dp), ALLOCATABLE, &
142 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 DIMENSION(:), POINTER :: ap_iterator
149 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
150 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
151 TYPE(gth_potential_type), POINTER :: gth_potential
152 REAL(kind=dp), DIMENSION(SIZE(particle_set)) :: at_thread
153 REAL(kind=dp), DIMENSION(nexp_max) :: alpha_ppl
154 REAL(kind=dp), DIMENSION(:, :), POINTER :: cval_lpot, h1_1block, h1_2block, &
155 h1_3block, h_block, p_block, rpgfa, &
156 rpgfb, sphi_a, sphi_b, zeta, zetb
157 REAL(kind=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
158 set_radius_a, set_radius_b
159 REAL(kind=dp), DIMENSION(4, nexp_max) :: cval_ppl
160 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!$ ALLOCATABLE, DIMENSION(:) :: locks
165!$ INTEGER :: lock_num, hash, hash1, hash2
166!$ INTEGER(KIND=int_8) :: iatom8
167!$ INTEGER, PARAMETER :: nlock = 501
168
169 do_dr = PRESENT(deltar)
170 doat = PRESENT(atcore)
171 IF ((calculate_forces .OR. doat) .AND. do_dr) THEN
172 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 libgrpp_local = .false.
179
180 IF (calculate_forces) THEN
181 CALL timeset(routinen//"_forces", handle)
182 ELSE
183 CALL timeset(routinen, handle)
184 END IF
185
186 nkind = SIZE(atomic_kind_set)
187 natom = SIZE(particle_set)
188
189 dokp = (nimages > 1)
190
191 IF (dokp) THEN
192 IF (PRESENT(cell_to_index)) THEN
193 cpassert(ASSOCIATED(cell_to_index))
194 ELSE
195 cpabort("Missing cell_to_index for k-point calculation")
196 END IF
197 END IF
198
199 IF (calculate_forces .OR. doat) THEN
200 IF (SIZE(matrix_p, 1) == 2) THEN
201 DO img = 1, nimages
202 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
203 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 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
206 END DO
207 END IF
208 END IF
209 force_thread = 0.0_dp
210 at_thread = 0.0_dp
211
212 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 basis_type=basis_type)
217
218 maxl = max(maxlgto, maxlppl)
219 CALL init_orbital_pointers(2*maxl + 2*nder + 1)
220
221 ldsab = max(maxco, ncoset(maxlppl), maxsgf, maxlppl)
222 ldai = ncoset(maxl + nder + 1)
223
224 ALLOCATE (basis_set_list(nkind))
225 DO ikind = 1, nkind
226 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
227 IF (ASSOCIATED(basis_set_a)) THEN
228 basis_set_list(ikind)%gto_basis_set => basis_set_a
229 ELSE
230 NULLIFY (basis_set_list(ikind)%gto_basis_set)
231 END IF
232 END DO
233
234 pv_thread = 0.0_dp
235
236 nthread = 1
237!$ nthread = omp_get_max_threads()
238
239 ! iterator for basis/potential list
240 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!$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 CALL neighbor_list_iterator_release(ap_iterator)
735
736 DEALLOCATE (basis_set_list)
737
738 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 IF (SIZE(matrix_p, 1) == 2) THEN
742 DO img = 1, nimages
743 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
744 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 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
747 END DO
748 END IF
749 END IF
750
751 IF (calculate_forces) THEN
752 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 atom_a = atom_of_kind(iatom)
756 ikind = kind_of(iatom)
757 force(ikind)%gth_ppl(:, atom_a) = force(ikind)%gth_ppl(:, atom_a) + force_thread(:, iatom)
758 END DO
759!$OMP END DO
760 DEALLOCATE (atom_of_kind, kind_of)
761 END IF
762 IF (doat) THEN
763 atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
764 END IF
765
766 IF (calculate_forces .AND. use_virial) THEN
767 virial%pv_ppl = virial%pv_ppl + pv_thread
768 virial%pv_virial = virial%pv_virial + pv_thread
769 END IF
770
771 CALL timestop(handle)
772
773 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
810 INTEGER, DIMENSION(1:10) :: nrloc
811 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, nct_lpot, npgfa, nsgfa
812 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: va, work
817 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: bcon, cval_lpot, rpgfa, sphi_a, zeta
826 REAL(kind=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
827 set_radius_a
828 REAL(kind=dp), DIMENSION(4, nexp_max) :: cval_ppl
829 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!$ ALLOCATABLE, DIMENSION(:) :: locks
834!$ INTEGER :: lock_num, hash
835!$ INTEGER, PARAMETER :: nlock = 501
836
837 IF (calculate_forces) THEN
838 CALL timeset(routinen//"_forces", handle)
839 ELSE
840 CALL timeset(routinen, handle)
841 END IF
842
843 nkind = SIZE(atomic_kind_set)
844 natom = SIZE(particle_set)
845
846 force_thread = 0.0_dp
847 pv_thread = 0.0_dp
848 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
849
850 ALLOCATE (basis_set_list(nkind))
851 DO ikind = 1, nkind
852 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
853 IF (ASSOCIATED(basis_set)) THEN
854 basis_set_list(ikind)%gto_basis_set => basis_set
855 ELSE
856 NULLIFY (basis_set_list(ikind)%gto_basis_set)
857 END IF
858 END DO
859
860 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!$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 IF (calculate_forces) THEN
1021 DO iatom = 1, natom
1022 atom_a = atom_of_kind(iatom)
1023 ikind = kind_of(iatom)
1024 force(ikind)%gth_ppl(1, atom_a) = force(ikind)%gth_ppl(1, atom_a) + force_thread(1, iatom)
1025 force(ikind)%gth_ppl(2, atom_a) = force(ikind)%gth_ppl(2, atom_a) + force_thread(2, iatom)
1026 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 DEALLOCATE (atom_of_kind, kind_of)
1030
1031 IF (calculate_forces .AND. use_virial) THEN
1032 virial%pv_ppl = virial%pv_ppl + pv_thread
1033 virial%pv_virial = virial%pv_virial + pv_thread
1034 END IF
1035
1036 DEALLOCATE (basis_set_list)
1037
1038 CALL timestop(handle)
1039
1040 END SUBROUTINE build_core_ppl_ri
1041
1042! **************************************************************************************************
1043
1044END MODULE core_ppl
external_potential_types::get_potential
Definition external_potential_types.F:278
qs_neighbor_list_types::nl_sub_iterate
Definition qs_neighbor_list_types.F:96
ai_overlap_ppl
Calculation of three-center overlap integrals over Cartesian Gaussian-type functions for the second t...
Definition ai_overlap_ppl.F:27
ai_overlap_ppl::ppl_integral
subroutine, public ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, hab2, hab2_work, deltar, iatom, jatom, katom)
Calculation of three-center overlap integrals <a|c|b> over Cartesian Gaussian functions for the local...
Definition ai_overlap_ppl.F:105
ai_overlap_ppl::ecploc_integral
subroutine, public ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, hab2, hab2_work, deltar, iatom, jatom, katom)
Calculation of three-center potential integrals <a|V(r)|b> over Cartesian Gaussian functions for the ...
Definition ai_overlap_ppl.F:239
ai_overlap_ppl::ppl_integral_ri
subroutine, public ppl_integral_ri(la_max_set, la_min_set, npgfa, rpgfa, zeta, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rac, dac, va, dva)
Calculation of two-center overlap integrals <a|c> over Cartesian Gaussian functions for the local par...
Definition ai_overlap_ppl.F:352
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
basis_set_types
Definition basis_set_types.F:15
core_ppl
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
Definition core_ppl.F:18
core_ppl::build_core_ppl_ri
subroutine, public build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_virial, qs_kind_set, atomic_kind_set, particle_set, sac_ppl, basis_type)
...
Definition core_ppl.F:791
core_ppl::build_core_ppl
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltar, atcore)
...
Definition core_ppl.F:97
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libgrpp_integrals
Local and semi-local ECP integrals using the libgrpp library.
Definition libgrpp_integrals.F:12
libgrpp_integrals::libgrpp_semilocal_integrals
subroutine, public libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Semi-local ECP integrals using libgrpp.
Definition libgrpp_integrals.F:275
libgrpp_integrals::libgrpp_local_integrals
subroutine, public libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Local ECP integrals using libgrpp.
Definition libgrpp_integrals.F:67
libgrpp_integrals::libgrpp_local_forces_ref
subroutine, public libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Reference local ECP force routine using l+-1 integrals. No call is made to the numerically unstable g...
Definition libgrpp_integrals.F:491
libgrpp_integrals::libgrpp_semilocal_forces_ref
subroutine, public libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
Reference semi-local ECP forces using l+-1 integrals. No call is made to the numerically unstable gra...
Definition libgrpp_integrals.F:660
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition lri_environment_types.F:18
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:486
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:979
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::nl_set_sub_iterator
subroutine, public nl_set_sub_iterator(iterator_set, ikind, jkind, iatom, mepos)
...
Definition qs_neighbor_list_types.F:288
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
external_potential_types::sgp_potential_type
Definition external_potential_types.F:174
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:332
particle_types::particle_type
Definition particle_types.F:35
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:185
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26