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