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