(git:b279b6b)
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 ! **************************************************************************************************
18 MODULE 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 
71 CONTAINS
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 
996 END MODULE core_ppl
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_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
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
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::nl_sub_iterate
integer function, public nl_sub_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:440
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