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