Line data Source code
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 17045 : 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 17045 : 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 17045 : 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 17045 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, &
126 17045 : nct_lpot, npgfa, npgfb, nsgfa, nsgfb
127 17045 : 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 17045 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qab, work
133 17045 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work
134 17045 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
135 : REAL(KIND=dp), ALLOCATABLE, &
136 17045 : 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 17045 : DIMENSION(:), POINTER :: ap_iterator
143 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
144 17045 : 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 17045 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h1_1block, h1_2block, &
148 17045 : h1_3block, h_block, p_block, rpgfa, &
149 17045 : rpgfb, sphi_a, sphi_b, zeta, zetb
150 17045 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
151 17045 : set_radius_a, set_radius_b
152 : REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
153 34090 : 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 17045 : !$ ALLOCATABLE, DIMENSION(:) :: locks
158 : !$ INTEGER :: lock_num, hash, hash1, hash2
159 : !$ INTEGER(KIND=int_8) :: iatom8
160 : !$ INTEGER, PARAMETER :: nlock = 501
161 :
162 17045 : do_dR = PRESENT(deltaR)
163 : MARK_USED(int_8)
164 :
165 34090 : IF (calculate_forces) THEN
166 7101 : CALL timeset(routineN//"_forces", handle)
167 : ELSE
168 9944 : CALL timeset(routineN, handle)
169 : END IF
170 :
171 17045 : nkind = SIZE(atomic_kind_set)
172 17045 : natom = SIZE(particle_set)
173 :
174 17045 : dokp = (nimages > 1)
175 :
176 17045 : IF (dokp) THEN
177 218 : CPASSERT(PRESENT(cell_to_index) .AND. ASSOCIATED(cell_to_index))
178 : END IF
179 :
180 17045 : IF (calculate_forces) THEN
181 7101 : IF (SIZE(matrix_p, 1) == 2) THEN
182 2390 : DO img = 1, nimages
183 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
184 1478 : 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 2390 : alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
187 : END DO
188 : END IF
189 : END IF
190 270837 : force_thread = 0.0_dp
191 :
192 17045 : 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 17045 : basis_type=basis_type)
197 :
198 17045 : maxl = MAX(maxlgto, maxlppl)
199 17045 : CALL init_orbital_pointers(2*maxl + 2*nder + 1)
200 :
201 17045 : ldsab = MAX(maxco, ncoset(maxlppl), maxsgf, maxlppl)
202 17045 : ldai = ncoset(maxl + nder + 1)
203 :
204 81302 : ALLOCATE (basis_set_list(nkind))
205 47212 : DO ikind = 1, nkind
206 30167 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
207 47212 : IF (ASSOCIATED(basis_set_a)) THEN
208 30167 : basis_set_list(ikind)%gto_basis_set => basis_set_a
209 : ELSE
210 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
211 : END IF
212 : END DO
213 :
214 17045 : pv_thread = 0.0_dp
215 :
216 : nthread = 1
217 17045 : !$ nthread = omp_get_max_threads()
218 :
219 : ! iterator for basis/potential list
220 17045 : 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 17045 : !$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 17045 : CALL neighbor_list_iterator_release(ap_iterator)
690 :
691 17045 : DEALLOCATE (basis_set_list)
692 :
693 17045 : 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 7101 : IF (SIZE(matrix_p, 1) == 2) THEN
697 2390 : DO img = 1, nimages
698 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
699 1478 : 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 2390 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
702 : END DO
703 : END IF
704 : END IF
705 :
706 17045 : IF (calculate_forces) THEN
707 7101 : 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 25121 : atom_a = atom_of_kind(iatom)
711 25121 : ikind = kind_of(iatom)
712 100484 : force(ikind)%gth_ppl(:, atom_a) = force(ikind)%gth_ppl(:, atom_a) + force_thread(:, iatom)
713 : END DO
714 : !$OMP END DO
715 7101 : DEALLOCATE (atom_of_kind, kind_of)
716 : END IF
717 :
718 17045 : IF (calculate_forces .AND. use_virial) THEN
719 10478 : virial%pv_ppl = virial%pv_ppl + pv_thread
720 10478 : virial%pv_virial = virial%pv_virial + pv_thread
721 : END IF
722 :
723 17045 : CALL timestop(handle)
724 :
725 51135 : 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 4 : 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 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
762 : INTEGER, DIMENSION(1:10) :: nrloc
763 4 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, nct_lpot, npgfa, nsgfa
764 4 : 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 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: va, work
769 4 : 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 4 : 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 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: bcon, cval_lpot, rpgfa, sphi_a, zeta
778 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
779 4 : set_radius_a
780 : REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
781 8 : 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 4 : !$ ALLOCATABLE, DIMENSION(:) :: locks
786 : !$ INTEGER :: lock_num, hash
787 : !$ INTEGER, PARAMETER :: nlock = 501
788 :
789 8 : IF (calculate_forces) THEN
790 2 : CALL timeset(routineN//"_forces", handle)
791 : ELSE
792 2 : CALL timeset(routineN, handle)
793 : END IF
794 :
795 4 : nkind = SIZE(atomic_kind_set)
796 4 : natom = SIZE(particle_set)
797 :
798 52 : force_thread = 0.0_dp
799 4 : pv_thread = 0.0_dp
800 4 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
801 :
802 20 : ALLOCATE (basis_set_list(nkind))
803 12 : DO ikind = 1, nkind
804 8 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
805 12 : IF (ASSOCIATED(basis_set)) THEN
806 8 : basis_set_list(ikind)%gto_basis_set => basis_set
807 : ELSE
808 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
809 : END IF
810 : END DO
811 :
812 4 : 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 4 : !$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 4 : IF (calculate_forces) THEN
973 8 : DO iatom = 1, natom
974 6 : atom_a = atom_of_kind(iatom)
975 6 : ikind = kind_of(iatom)
976 6 : force(ikind)%gth_ppl(1, atom_a) = force(ikind)%gth_ppl(1, atom_a) + force_thread(1, iatom)
977 6 : force(ikind)%gth_ppl(2, atom_a) = force(ikind)%gth_ppl(2, atom_a) + force_thread(2, iatom)
978 8 : 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 4 : DEALLOCATE (atom_of_kind, kind_of)
982 :
983 4 : IF (calculate_forces .AND. use_virial) THEN
984 0 : virial%pv_ppl = virial%pv_ppl + pv_thread
985 0 : virial%pv_virial = virial%pv_virial + pv_thread
986 : END IF
987 :
988 4 : DEALLOCATE (basis_set_list)
989 :
990 4 : CALL timestop(handle)
991 :
992 12 : END SUBROUTINE build_core_ppl_ri
993 :
994 : ! **************************************************************************************************
995 :
996 : END MODULE core_ppl
|
