Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief Calculation of the local pseudopotential contribution to the core Hamiltonian
9 : !> <a|V(local)|b> = <a|Sum e^a*rc**2|b>
10 : !> \par History
11 : !> - core_ppnl refactored from qs_core_hamiltonian [Joost VandeVondele, 2008-11-01]
12 : !> - adapted for PPL [jhu, 2009-02-23]
13 : !> - OpenMP added [Iain Bethune, Fiona Reid, 2013-11-13]
14 : !> - Bug fix: correct orbital pointer range [07.2014,JGH]
15 : !> - k-point aware [07.2015,JGH]
16 : !> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
17 : ! **************************************************************************************************
18 : MODULE core_ppl
19 :
20 : USE ai_overlap_ppl, ONLY: ecploc_integral,&
21 : ppl_integral,&
22 : ppl_integral_ri
23 : USE atomic_kind_types, ONLY: atomic_kind_type,&
24 : get_atomic_kind_set
25 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
26 : gto_basis_set_type
27 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
28 : dbcsr_get_block_p,&
29 : dbcsr_p_type
30 : USE external_potential_types, ONLY: get_potential,&
31 : gth_potential_type,&
32 : sgp_potential_type
33 : USE kinds, ONLY: dp,&
34 : int_8
35 : USE libgrpp_integrals, ONLY: libgrpp_local_forces_ref,&
36 : libgrpp_local_integrals,&
37 : libgrpp_semilocal_forces_ref,&
38 : libgrpp_semilocal_integrals
39 : USE lri_environment_types, ONLY: lri_kind_type
40 : USE orbital_pointers, ONLY: init_orbital_pointers,&
41 : ncoset
42 : USE particle_types, ONLY: particle_type
43 : USE qs_force_types, ONLY: qs_force_type
44 : USE qs_kind_types, ONLY: get_qs_kind,&
45 : get_qs_kind_set,&
46 : qs_kind_type
47 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
48 : neighbor_list_iterator_create,&
49 : neighbor_list_iterator_p_type,&
50 : neighbor_list_iterator_release,&
51 : neighbor_list_set_p_type,&
52 : nl_set_sub_iterator,&
53 : nl_sub_iterate
54 : USE upf_local_integrals, ONLY: upf_local_init,&
55 : upf_local_matrix,&
56 : upf_local_ri
57 : USE virial_methods, ONLY: virial_pair_force
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 :
73 : PUBLIC :: build_core_ppl, build_core_ppl_ri
74 :
75 : CONTAINS
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 21070 : 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 21070 : 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 21070 : 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 21070 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, &
133 21070 : nct_lpot, npgfa, npgfb, nsgfa, nsgfb
134 21070 : 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 21070 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: upf_v, work
143 21070 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work, upf_da, &
144 21070 : upf_db
145 21070 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
146 : REAL(KIND=dp), ALLOCATABLE, &
147 21070 : 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
152 : TYPE(neighbor_list_iterator_p_type), &
153 21070 : DIMENSION(:), POINTER :: ap_iterator
154 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
155 21070 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
156 : TYPE(gth_potential_type), POINTER :: gth_potential
157 42140 : REAL(KIND=dp), DIMENSION(SIZE(particle_set)) :: at_thread
158 : REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
159 21070 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h1_1block, h1_2block, &
160 21070 : h1_3block, h_block, p_block, rpgfa, &
161 21070 : rpgfb, sphi_a, sphi_b, zeta, zetb
162 21070 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
163 21070 : set_radius_a, set_radius_b
164 : REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
165 42140 : 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 21070 : !$ ALLOCATABLE, DIMENSION(:) :: locks
170 : !$ INTEGER :: lock_num, hash, hash1, hash2
171 : !$ INTEGER(KIND=int_8) :: iatom8
172 : !$ INTEGER, PARAMETER :: nlock = 501
173 :
174 21070 : do_dR = PRESENT(deltaR)
175 21070 : doat = PRESENT(atcore)
176 21070 : IF ((calculate_forces .OR. doat) .AND. do_dR) THEN
177 0 : 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 21070 : libgrpp_local = .FALSE.
184 :
185 21070 : IF (calculate_forces) THEN
186 7977 : CALL timeset(routineN//"_forces", handle)
187 : ELSE
188 13093 : CALL timeset(routineN, handle)
189 : END IF
190 :
191 21070 : nkind = SIZE(atomic_kind_set)
192 21070 : natom = SIZE(particle_set)
193 :
194 21070 : dokp = (nimages > 1)
195 :
196 21070 : IF (dokp) THEN
197 1272 : IF (PRESENT(cell_to_index)) THEN
198 1272 : CPASSERT(ASSOCIATED(cell_to_index))
199 : ELSE
200 0 : CPABORT("Missing cell_to_index for k-point calculation")
201 : END IF
202 : END IF
203 :
204 21070 : IF (calculate_forces .OR. doat) THEN
205 8039 : IF (SIZE(matrix_p, 1) == 2) THEN
206 3734 : DO img = 1, nimages
207 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
208 2750 : 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 3734 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
211 : END DO
212 : END IF
213 : END IF
214 331030 : force_thread = 0.0_dp
215 98560 : at_thread = 0.0_dp
216 :
217 21070 : 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 21070 : basis_type=basis_type)
222 :
223 21070 : maxl = MAX(maxlgto, maxlppl)
224 21070 : CALL init_orbital_pointers(2*maxl + 2*nder + 1)
225 :
226 21070 : ldsab = MAX(maxco, ncoset(maxlppl), maxsgf, maxlppl)
227 21070 : ldai = ncoset(maxl + nder + 1)
228 :
229 100011 : ALLOCATE (basis_set_list(nkind))
230 57871 : DO ikind = 1, nkind
231 36801 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
232 57871 : IF (ASSOCIATED(basis_set_a)) THEN
233 36801 : basis_set_list(ikind)%gto_basis_set => basis_set_a
234 : ELSE
235 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
236 : END IF
237 : END DO
238 :
239 21070 : CALL prepare_upf_local(qs_kind_set, 2*maxlgto + 1, have_upf)
240 :
241 21070 : pv_thread = 0.0_dp
242 :
243 : nthread = 1
244 21070 : !$ nthread = omp_get_max_threads()
245 :
246 : ! iterator for basis/potential list
247 21070 : 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 21070 : !$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 21070 : CALL neighbor_list_iterator_release(ap_iterator)
781 :
782 21070 : DEALLOCATE (basis_set_list)
783 :
784 21070 : 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 8039 : IF (SIZE(matrix_p, 1) == 2) THEN
788 3734 : DO img = 1, nimages
789 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
790 2750 : 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 3734 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
793 : END DO
794 : END IF
795 : END IF
796 :
797 21070 : IF (calculate_forces) THEN
798 7977 : 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 28335 : atom_a = atom_of_kind(iatom)
802 28335 : ikind = kind_of(iatom)
803 113340 : force(ikind)%gth_ppl(:, atom_a) = force(ikind)%gth_ppl(:, atom_a) + force_thread(:, iatom)
804 : END DO
805 : !$OMP END DO
806 7977 : DEALLOCATE (atom_of_kind, kind_of)
807 : END IF
808 21070 : IF (doat) THEN
809 280 : atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
810 : END IF
811 :
812 21070 : IF (calculate_forces .AND. use_virial) THEN
813 12558 : virial%pv_ppl = virial%pv_ppl + pv_thread
814 12558 : virial%pv_virial = virial%pv_virial + pv_thread
815 : END IF
816 :
817 21070 : CALL timestop(handle)
818 :
819 63210 : 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 4 : 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 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
856 : INTEGER, DIMENSION(1:10) :: nrloc
857 4 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, nct_lpot, npgfa, nsgfa
858 4 : 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 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: va, work
864 4 : 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 4 : 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 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: bcon, cval_lpot, rpgfa, sphi_a, zeta
873 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
874 4 : set_radius_a
875 : REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
876 8 : 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 4 : !$ ALLOCATABLE, DIMENSION(:) :: locks
881 : !$ INTEGER :: lock_num, hash
882 : !$ INTEGER, PARAMETER :: nlock = 501
883 :
884 8 : IF (calculate_forces) THEN
885 2 : CALL timeset(routineN//"_forces", handle)
886 : ELSE
887 2 : CALL timeset(routineN, handle)
888 : END IF
889 :
890 4 : nkind = SIZE(atomic_kind_set)
891 4 : natom = SIZE(particle_set)
892 :
893 52 : force_thread = 0.0_dp
894 4 : pv_thread = 0.0_dp
895 4 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
896 :
897 20 : ALLOCATE (basis_set_list(nkind))
898 12 : DO ikind = 1, nkind
899 8 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
900 12 : IF (ASSOCIATED(basis_set)) THEN
901 8 : basis_set_list(ikind)%gto_basis_set => basis_set
902 : ELSE
903 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
904 : END IF
905 : END DO
906 :
907 4 : CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxsgf=maxsgf, maxlgto=maxlgto, basis_type=basis_type)
908 4 : 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 4 : !$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 4 : IF (calculate_forces) THEN
1081 8 : DO iatom = 1, natom
1082 6 : atom_a = atom_of_kind(iatom)
1083 6 : ikind = kind_of(iatom)
1084 6 : force(ikind)%gth_ppl(1, atom_a) = force(ikind)%gth_ppl(1, atom_a) + force_thread(1, iatom)
1085 6 : force(ikind)%gth_ppl(2, atom_a) = force(ikind)%gth_ppl(2, atom_a) + force_thread(2, iatom)
1086 8 : 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 4 : DEALLOCATE (atom_of_kind, kind_of)
1090 :
1091 4 : IF (calculate_forces .AND. use_virial) THEN
1092 0 : virial%pv_ppl = virial%pv_ppl + pv_thread
1093 0 : virial%pv_virial = virial%pv_virial + pv_thread
1094 : END IF
1095 :
1096 4 : DEALLOCATE (basis_set_list)
1097 :
1098 4 : CALL timestop(handle)
1099 :
1100 12 : 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 21074 : 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 21074 : IF (PRESENT(have_upf)) have_upf = .FALSE.
1119 57883 : DO ikind = 1, SIZE(qs_kind_set)
1120 36809 : CALL get_qs_kind(qs_kind_set(ikind), sgp_potential=potential)
1121 36809 : IF (.NOT. ASSOCIATED(potential)) CYCLE
1122 234 : IF (.NOT. ALLOCATED(potential%upf_data)) CYCLE
1123 36 : IF (PRESENT(have_upf)) have_upf = .TRUE.
1124 21110 : IF (potential%upf_local%radial%degree < degree) THEN
1125 24 : CALL upf_local_init(potential%upf_data, potential%alpha_core_charge, degree, potential%upf_local)
1126 : END IF
1127 : END DO
1128 21074 : END SUBROUTINE prepare_upf_local
1129 :
1130 : END MODULE core_ppl
|
