(git:3f0d10b)
Loading...
Searching...
No Matches
core_ae.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Calculation of the nuclear attraction contribution to the core Hamiltonian
9!> <a|erfc|b> :we only calculate the non-screened part
10!> \par History
11!> - core_ppnl refactored from qs_core_hamiltonian [Joost VandeVondele, 2008-11-01]
12!> - adapted for nuclear attraction [jhu, 2009-02-24]
13! **************************************************************************************************
14MODULE core_ae
15 USE ai_verfc, ONLY: verfc
16 USE ao_util, ONLY: exp_radius
17 USE atomic_kind_types, ONLY: atomic_kind_type,&
18 get_atomic_kind_set
19 USE basis_set_types, ONLY: gto_basis_set_p_type,&
20 gto_basis_set_type
21 USE cp_dbcsr_api, ONLY: dbcsr_add,&
22 dbcsr_get_block_p,&
23 dbcsr_p_type
24 USE external_potential_types, ONLY: all_potential_type,&
25 get_potential,&
26 sgp_potential_type
27 USE kinds, ONLY: dp,&
28 int_8
29 USE orbital_pointers, ONLY: coset,&
30 indco,&
31 init_orbital_pointers,&
32 ncoset
33 USE particle_types, ONLY: particle_type
34 USE qs_force_types, ONLY: qs_force_type
35 USE qs_kind_types, ONLY: get_qs_kind,&
36 get_qs_kind_set,&
37 qs_kind_type
38 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
39 neighbor_list_iterator_create,&
40 neighbor_list_iterator_p_type,&
41 neighbor_list_iterator_release,&
42 neighbor_list_set_p_type,&
43 nl_set_sub_iterator,&
44 nl_sub_iterate
45 USE virial_methods, ONLY: virial_pair_force
46 USE virial_types, ONLY: virial_type
47
48!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
49!$ USE OMP_LIB, ONLY: omp_lock_kind, &
50!$ omp_init_lock, omp_set_lock, &
51!$ omp_unset_lock, omp_destroy_lock
52
53#include "./base/base_uses.f90"
54
55 IMPLICIT NONE
56
57 PRIVATE
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ae'
60
61 PUBLIC :: build_core_ae, build_erfc, verfc_force
62
63CONTAINS
64
65! **************************************************************************************************
66!> \brief ...
67!> \param matrix_h ...
68!> \param matrix_p ...
69!> \param force ...
70!> \param virial ...
71!> \param calculate_forces ...
72!> \param use_virial ...
73!> \param nder ...
74!> \param qs_kind_set ...
75!> \param atomic_kind_set ...
76!> \param particle_set ...
77!> \param sab_orb ...
78!> \param sac_ae ...
79!> \param nimages ...
80!> \param cell_to_index ...
81!> \param basis_type ...
82!> \param atcore ...
83! **************************************************************************************************
84 SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
85 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
86 nimages, cell_to_index, basis_type, atcore)
87
88 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
89 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
90 TYPE(virial_type), POINTER :: virial
91 LOGICAL, INTENT(IN) :: calculate_forces
92 LOGICAL :: use_virial
93 INTEGER :: nder
94 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
95 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
96 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
97 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
98 POINTER :: sab_orb, sac_ae
99 INTEGER, INTENT(IN) :: nimages
100 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
101 CHARACTER(LEN=*), INTENT(IN) :: basis_type
102 REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
103 OPTIONAL :: atcore
104
105 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_core_ae'
106
107 INTEGER :: atom_a, handle, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, katom, &
108 kkind, ldai, ldsab, max_na_plus, maxco, maxl, maxnset, maxsgf, mepos, na_plus, natom, &
109 nb_plus, ncoa, ncob, nij, nkind, nseta, nsetb, nthread, sgfa, sgfb, slot
110 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
111 INTEGER, DIMENSION(3) :: cellind
112 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
113 npgfb, nsgfa, nsgfb
114 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
115 LOGICAL :: any_habd, doat, dokp, found, need_pab
116 LOGICAL, ALLOCATABLE, DIMENSION(:) :: active_set_pair, has_core_potential, &
117 pab_ready
118 REAL(kind=dp) :: alpha_c, atk0, atk1, core_charge, core_radius, dab, dac, dbc, f0, rab2, &
119 rac2, rbc2, set_radius_a_max, set_radius_b_max, zeta_c
120 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: alpha_core_kind, core_charge_kind, &
121 core_radius_kind, ff, zeta_core_kind
122 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
123 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
124 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
125 REAL(kind=dp), DIMENSION(3, 3) :: pv_thread
126 TYPE(neighbor_list_iterator_p_type), &
127 DIMENSION(:), POINTER :: ap_iterator
128 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
129 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
130 TYPE(all_potential_type), POINTER :: all_potential
131 REAL(kind=dp), DIMENSION(SIZE(particle_set)) :: at_thread
132 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
133 sphi_b, zeta, zetb
134 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
135 REAL(kind=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
136 TYPE(sgp_potential_type), POINTER :: sgp_potential
137
138!$ INTEGER(kind=omp_lock_kind), &
139!$ ALLOCATABLE, DIMENSION(:) :: locks
140!$ INTEGER :: lock_num, hash, hash1, hash2
141!$ INTEGER(KIND=int_8) :: iatom8
142!$ INTEGER, PARAMETER :: nlock = 501
143
144 mark_used(int_8)
145
146 IF (calculate_forces) THEN
147 CALL timeset(routinen//"_forces", handle)
148 ELSE
149 CALL timeset(routinen, handle)
150 END IF
151
152 nkind = SIZE(atomic_kind_set)
153 natom = SIZE(particle_set)
154
155 doat = PRESENT(atcore)
156 dokp = (nimages > 1)
157
158 IF (calculate_forces .OR. doat) THEN
159 IF (SIZE(matrix_p, 1) == 2) THEN
160 DO img = 1, nimages
161 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
162 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
163 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
164 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
165 END DO
166 END IF
167 END IF
168
169 force_thread = 0.0_dp
170 at_thread = 0.0_dp
171 pv_thread = 0.0_dp
172
173 ALLOCATE (basis_set_list(nkind))
174 ALLOCATE (has_core_potential(nkind), alpha_core_kind(nkind), zeta_core_kind(nkind), &
175 core_charge_kind(nkind), core_radius_kind(nkind))
176 has_core_potential(:) = .false.
177 alpha_core_kind(:) = 0.0_dp
178 zeta_core_kind(:) = 0.0_dp
179 core_charge_kind(:) = 0.0_dp
180 core_radius_kind(:) = 0.0_dp
181 DO ikind = 1, nkind
182 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type)
183 IF (ASSOCIATED(basis_set_a)) THEN
184 basis_set_list(ikind)%gto_basis_set => basis_set_a
185 ELSE
186 NULLIFY (basis_set_list(ikind)%gto_basis_set)
187 END IF
188 CALL get_qs_kind(qs_kind_set(ikind), all_potential=all_potential, &
189 sgp_potential=sgp_potential)
190 IF (ASSOCIATED(all_potential)) THEN
191 CALL get_potential(potential=all_potential, &
192 alpha_core_charge=alpha_core_kind(ikind), zeff=zeta_core_kind(ikind), &
193 ccore_charge=core_charge_kind(ikind), core_charge_radius=core_radius_kind(ikind))
194 has_core_potential(ikind) = .true.
195 ELSE IF (ASSOCIATED(sgp_potential)) THEN
196 CALL get_potential(potential=sgp_potential, &
197 alpha_core_charge=alpha_core_kind(ikind), zeff=zeta_core_kind(ikind), &
198 ccore_charge=core_charge_kind(ikind), core_charge_radius=core_radius_kind(ikind))
199 has_core_potential(ikind) = .true.
200 END IF
201 END DO
202
203 CALL get_qs_kind_set(qs_kind_set, basis_type=basis_type, &
204 maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset)
205 max_na_plus = 0
206 IF (calculate_forces) CALL get_qs_kind_set(qs_kind_set, basis_type=basis_type, &
207 maxco=max_na_plus, maxder=nder)
208 CALL init_orbital_pointers(maxl + nder + 1)
209 ldsab = max(maxco, maxsgf)
210 ldai = ncoset(maxl + nder + 1)
211
212 nthread = 1
213!$ nthread = omp_get_max_threads()
214
215 ! iterator for basis/potential list
216 CALL neighbor_list_iterator_create(ap_iterator, sac_ae, search=.true., nthread=nthread)
217
218!$OMP PARALLEL &
219!$OMP DEFAULT (NONE) &
220!$OMP SHARED (ap_iterator, basis_set_list, calculate_forces, use_virial, &
221!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
222!$OMP sab_orb, sac_ae, nthread, ncoset, nkind, cell_to_index, &
223!$OMP slot, ldsab, maxnset, ldai, nder, maxl, maxco, max_na_plus, dokp, doat, locks, &
224!$OMP natom, has_core_potential, alpha_core_kind, zeta_core_kind, core_charge_kind, &
225!$OMP core_radius_kind) &
226!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
227!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, lock_num, &
228!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
229!$OMP zetb, zeta_c, alpha_c, core_charge, dab, irow, icol, h_block, found, iset, ncoa, &
230!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
231!$OMP rac, dac, rbc, rab2, rac2, rbc2, dbc, na_plus, nb_plus, verf, vnuc, &
232!$OMP set_radius_a, core_radius, set_radius_a_max, set_radius_b_max, rpgfa, force_a, force_b, mepos, &
233!$OMP atk0, atk1, habd, f0, katom, cellind, img, nij, ff, any_habd, need_pab, pab_ready, &
234!$OMP active_set_pair, sgp_potential, all_potential, hash, hash1, hash2, iatom8) &
235!$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )
236
237!$OMP SINGLE
238!$ ALLOCATE (locks(nlock))
239!$OMP END SINGLE
240
241!$OMP DO
242!$ DO lock_num = 1, nlock
243!$ call omp_init_lock(locks(lock_num))
244!$ END DO
245!$OMP END DO
246
247 mepos = 0
248!$ mepos = omp_get_thread_num()
249
250 ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab))
251 ALLOCATE (verf(ldai, ldai, 2*maxl + nder + 1), vnuc(ldai, ldai, 2*maxl + nder + 1), ff(0:2*maxl + nder))
252 ALLOCATE (active_set_pair(maxnset*maxnset))
253 IF (calculate_forces .OR. doat) ALLOCATE (pab(maxco, maxco, maxnset*maxnset), pab_ready(maxnset*maxnset))
254 IF (calculate_forces) ALLOCATE (habd(max(1, max_na_plus), max(1, maxco)))
255
256!$OMP DO SCHEDULE(GUIDED)
257 DO slot = 1, sab_orb(1)%nl_size
258
259 ikind = sab_orb(1)%nlist_task(slot)%ikind
260 jkind = sab_orb(1)%nlist_task(slot)%jkind
261 iatom = sab_orb(1)%nlist_task(slot)%iatom
262 jatom = sab_orb(1)%nlist_task(slot)%jatom
263 cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
264 rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
265
266 basis_set_a => basis_set_list(ikind)%gto_basis_set
267 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
268 basis_set_b => basis_set_list(jkind)%gto_basis_set
269 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
270!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
271!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
272 ! basis ikind
273 first_sgfa => basis_set_a%first_sgf
274 la_max => basis_set_a%lmax
275 la_min => basis_set_a%lmin
276 npgfa => basis_set_a%npgf
277 nseta = basis_set_a%nset
278 nsgfa => basis_set_a%nsgf_set
279 rpgfa => basis_set_a%pgf_radius
280 set_radius_a => basis_set_a%set_radius
281 sphi_a => basis_set_a%sphi
282 zeta => basis_set_a%zet
283 set_radius_a_max = maxval(set_radius_a(1:nseta))
284 ! basis jkind
285 first_sgfb => basis_set_b%first_sgf
286 lb_max => basis_set_b%lmax
287 lb_min => basis_set_b%lmin
288 npgfb => basis_set_b%npgf
289 nsetb = basis_set_b%nset
290 nsgfb => basis_set_b%nsgf_set
291 rpgfb => basis_set_b%pgf_radius
292 set_radius_b => basis_set_b%set_radius
293 sphi_b => basis_set_b%sphi
294 zetb => basis_set_b%zet
295 set_radius_b_max = maxval(set_radius_b(1:nsetb))
296
297 dab = sqrt(sum(rab*rab))
298 IF (set_radius_a_max + set_radius_b_max < dab) cycle
299 rab2 = dab*dab
300
301 DO iset = 1, nseta
302 DO jset = 1, nsetb
303 nij = jset + (iset - 1)*maxnset
304 active_set_pair(nij) = set_radius_a(iset) + set_radius_b(jset) >= dab
305 END DO
306 END DO
307
308 IF (dokp) THEN
309 img = cell_to_index(cellind(1), cellind(2), cellind(3))
310 ELSE
311 img = 1
312 END IF
313
314 ! *** Use the symmetry of the first derivatives ***
315 IF (iatom == jatom) THEN
316 f0 = 1.0_dp
317 ELSE
318 f0 = 2.0_dp
319 END IF
320
321 ! *** Create matrix blocks for a new matrix block column ***
322 IF (iatom <= jatom) THEN
323 irow = iatom
324 icol = jatom
325 ELSE
326 irow = jatom
327 icol = iatom
328 END IF
329 NULLIFY (h_block)
330 CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, &
331 row=irow, col=icol, block=h_block, found=found)
332 IF (calculate_forces .OR. doat) THEN
333 NULLIFY (p_block)
334 CALL dbcsr_get_block_p(matrix=matrix_p(1, img)%matrix, &
335 row=irow, col=icol, block=p_block, found=found)
336 cpassert(ASSOCIATED(p_block))
337 pab_ready(1:maxnset*maxnset) = .false.
338 END IF
339
340 ! loop over all kinds for pseudopotential atoms
341 DO iset = 1, nseta
342 ncoa = npgfa(iset)*ncoset(la_max(iset))
343 DO jset = 1, nsetb
344 nij = jset + (iset - 1)*maxnset
345 IF (.NOT. active_set_pair(nij)) cycle
346 ncob = npgfb(jset)*ncoset(lb_max(jset))
347 hab(1:ncoa, 1:ncob, nij) = 0._dp
348 END DO
349 END DO
350 DO kkind = 1, nkind
351 IF (.NOT. has_core_potential(kkind)) cycle
352 alpha_c = alpha_core_kind(kkind)
353 zeta_c = zeta_core_kind(kkind)
354 core_charge = core_charge_kind(kkind)
355 core_radius = core_radius_kind(kkind)
356
357 CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
358
359 DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
360 CALL get_iterator_info(ap_iterator, jatom=katom, r=rac, mepos=mepos)
361
362 dac = sqrt(sum(rac*rac))
363 rbc(:) = rac(:) - rab(:)
364 dbc = sqrt(sum(rbc*rbc))
365 rac2 = dac*dac
366 rbc2 = dbc*dbc
367 IF ((set_radius_a_max + core_radius < dac) .OR. &
368 (set_radius_b_max + core_radius < dbc)) THEN
369 cycle
370 END IF
371
372 DO iset = 1, nseta
373 IF (set_radius_a(iset) + core_radius < dac) cycle
374 ncoa = npgfa(iset)*ncoset(la_max(iset))
375 sgfa = first_sgfa(1, iset)
376 DO jset = 1, nsetb
377 nij = jset + (iset - 1)*maxnset
378 IF (.NOT. active_set_pair(nij)) cycle
379 IF (set_radius_b(jset) + core_radius < dbc) cycle
380 ncob = npgfb(jset)*ncoset(lb_max(jset))
381 sgfb = first_sgfb(1, jset)
382 ! *** Calculate the GTH pseudo potential forces ***
383 atk0 = 0.0_dp
384 IF (doat) THEN
385 IF (pab_ready(nij)) atk0 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
386 END IF
387 IF (calculate_forces) THEN
388 na_plus = npgfa(iset)*ncoset(la_max(iset) + nder)
389 nb_plus = npgfb(jset)*ncoset(lb_max(jset))
390 habd(1:na_plus, 1:nb_plus) = 0._dp
391 CALL verfc( &
392 la_max(iset) + nder, npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
393 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
394 alpha_c, core_radius, zeta_c, core_charge, &
395 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:), &
396 nder, habd)
397 any_habd = any(habd(1:na_plus, 1:nb_plus) /= 0.0_dp)
398 ELSE
399 CALL verfc( &
400 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
401 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
402 alpha_c, core_radius, zeta_c, core_charge, &
403 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:))
404 any_habd = .false.
405 END IF
406 IF (calculate_forces .OR. doat) THEN
407 IF (.NOT. pab_ready(nij)) THEN
408 need_pab = any_habd
409 IF (doat .AND. .NOT. need_pab) need_pab = any(hab(1:ncoa, 1:ncob, nij) /= 0.0_dp)
410 IF (need_pab) THEN
411 ! *** Decontract density matrix block ***
412 IF (iatom <= jatom) THEN
413 work(1:ncoa, 1:nsgfb(jset)) = &
414 matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
415 p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
416 ELSE
417 work(1:ncoa, 1:nsgfb(jset)) = &
418 matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
419 transpose(p_block(sgfb:sgfb + nsgfb(jset) - 1, &
420 sgfa:sgfa + nsgfa(iset) - 1)))
421 END IF
422 pab(1:ncoa, 1:ncob, nij) = &
423 matmul(work(1:ncoa, 1:nsgfb(jset)), &
424 transpose(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
425 pab_ready(nij) = .true.
426 END IF
427 END IF
428 END IF
429 IF (any_habd) THEN
430 ! *** The derivatives w.r.t. atomic center c are ***
431 ! *** calculated using the translational invariance ***
432 ! *** of the first derivatives ***
433 CALL verfc_force(habd, pab(:, :, nij), force_a, force_b, nder, &
434 la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
435 lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), rab)
436 force_thread(1:3, iatom) = force_thread(1:3, iatom) + f0*force_a(1:3)
437 force_thread(1:3, jatom) = force_thread(1:3, jatom) + f0*force_b(1:3)
438 force_thread(1:3, katom) = force_thread(1:3, katom) - f0*force_a(1:3) - f0*force_b(1:3)
439 IF (use_virial) THEN
440 CALL virial_pair_force(pv_thread, f0, force_a, rac)
441 CALL virial_pair_force(pv_thread, f0, force_b, rbc)
442 END IF
443 END IF
444 ! calculate atomic contributions
445 IF (doat) THEN
446 IF (pab_ready(nij)) THEN
447 atk1 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
448 at_thread(katom) = at_thread(katom) + (atk1 - atk0)
449 END IF
450 END IF
451 END DO
452 END DO
453 END DO
454 END DO
455 ! *** Contract nuclear attraction integrals
456 DO iset = 1, nseta
457 sgfa = first_sgfa(1, iset)
458 ncoa = npgfa(iset)*ncoset(la_max(iset))
459 DO jset = 1, nsetb
460 nij = jset + (iset - 1)*maxnset
461 IF (.NOT. active_set_pair(nij)) cycle
462 ncob = npgfb(jset)*ncoset(lb_max(jset))
463 IF (.NOT. any(hab(1:ncoa, 1:ncob, nij) /= 0.0_dp)) cycle
464 sgfb = first_sgfb(1, jset)
465!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
466!$ hash = MOD(hash1 + hash2, nlock) + 1
467 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab(1:ncoa, 1:ncob, nij), &
468 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
469!$ CALL omp_set_lock(locks(hash))
470 IF (iatom <= jatom) THEN
471 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
472 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
473 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
474 ELSE
475 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
476 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
477 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
478 END IF
479!$ CALL omp_unset_lock(locks(hash))
480 END DO
481 END DO
482
483 END DO
484
485 DEALLOCATE (hab, work, verf, vnuc, ff, active_set_pair)
486 IF (calculate_forces .OR. doat) DEALLOCATE (pab, pab_ready)
487 IF (calculate_forces) DEALLOCATE (habd)
488
489!$OMP DO
490!$ DO lock_num = 1, nlock
491!$ call omp_destroy_lock(locks(lock_num))
492!$ END DO
493!$OMP END DO
494
495!$OMP SINGLE
496!$ DEALLOCATE (locks)
497!$OMP END SINGLE NOWAIT
498
499!$OMP END PARALLEL
500
501 CALL neighbor_list_iterator_release(ap_iterator)
502
503 DEALLOCATE (basis_set_list, has_core_potential, alpha_core_kind, zeta_core_kind, &
504 core_charge_kind, core_radius_kind)
505
506 IF (calculate_forces .OR. doat) THEN
507 ! *** If LSD, then recover alpha density and beta density ***
508 ! *** from the total density (1) and the spin density (2) ***
509 IF (SIZE(matrix_p, 1) == 2) THEN
510 DO img = 1, nimages
511 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
512 alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
513 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
514 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
515 END DO
516 END IF
517 END IF
518
519 IF (calculate_forces) THEN
520 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
521!$OMP DO
522 DO iatom = 1, natom
523 atom_a = atom_of_kind(iatom)
524 ikind = kind_of(iatom)
525 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) + force_thread(:, iatom)
526 END DO
527!$OMP END DO
528 END IF
529 IF (doat) THEN
530 atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
531 END IF
532
533 IF (calculate_forces .AND. use_virial) THEN
534 virial%pv_ppl = virial%pv_ppl + pv_thread
535 virial%pv_virial = virial%pv_virial + pv_thread
536 END IF
537
538 CALL timestop(handle)
539
540 END SUBROUTINE build_core_ae
541
542! **************************************************************************************************
543!> \brief Forces from [a|V_erfc|b] nuclear attraction integrals for one shell pair
544!> \param habd Primitive [a|V_erfc|b] integrals up to la_max + nder (stride ncoset(la_max + nder))
545!> \param pab Decontracted density matrix block for the shell pair
546!> \param fa Force contribution on center A
547!> \param fb Force contribution on center B
548!> \param nder Derivative order (nder >= 1)
549!> \param la_max Maximum angular momentum on center A
550!> \param la_min Minimum angular momentum on center A
551!> \param npgfa Number of primitives on center A
552!> \param zeta Exponents on center A
553!> \param lb_max Maximum angular momentum on center B
554!> \param lb_min Minimum angular momentum on center B
555!> \param npgfb Number of primitives on center B
556!> \param zetb Exponents on center B
557!> \param rab Distance vector between centers A and B (R_B - R_A)
558! **************************************************************************************************
559 SUBROUTINE verfc_force(habd, pab, fa, fb, nder, la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, rab)
560
561 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: habd, pab
562 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: fa, fb
563 INTEGER, INTENT(IN) :: nder, la_max, la_min, npgfa
564 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
565 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
566 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
567 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
568
569 INTEGER :: ic_a, ic_a_min, ic_b, ic_b_min, icax, &
570 icbx, icob, icob_m, icpa, idir, ipgfa, &
571 ipgfb, na, nap, nb
572 INTEGER, DIMENSION(3) :: la, lb
573 INTEGER, DIMENSION(ncoset(lb_max), 3) :: icbm
574 INTEGER, DIMENSION(ncoset(la_max), 3) :: icam, icap
575 REAL(kind=dp) :: pval, s0, smb, zax2, zbx2
576 REAL(kind=dp), DIMENSION(3) :: sma, sp
577 REAL(kind=dp), DIMENSION(ncoset(la_max), 3) :: rla
578
579 fa = 0.0_dp
580 fb = 0.0_dp
581
582 na = ncoset(la_max)
583 nap = ncoset(la_max + nder)
584 nb = ncoset(lb_max)
585 ic_a_min = ncoset(la_min - 1) + 1
586 ic_b_min = ncoset(lb_min - 1) + 1
587 DO ic_a = ic_a_min, na
588 la(1:3) = indco(1:3, ic_a)
589 rla(ic_a, 1:3) = real(la(1:3), kind=dp)
590 icap(ic_a, 1) = coset(la(1) + 1, la(2), la(3))
591 icap(ic_a, 2) = coset(la(1), la(2) + 1, la(3))
592 icap(ic_a, 3) = coset(la(1), la(2), la(3) + 1)
593 icam(ic_a, 1) = coset(la(1) - 1, la(2), la(3))
594 icam(ic_a, 2) = coset(la(1), la(2) - 1, la(3))
595 icam(ic_a, 3) = coset(la(1), la(2), la(3) - 1)
596 END DO
597 DO ic_b = ic_b_min, nb
598 lb(1:3) = indco(1:3, ic_b)
599 icbm(ic_b, 1) = coset(lb(1) - 1, lb(2), lb(3))
600 icbm(ic_b, 2) = coset(lb(1), lb(2) - 1, lb(3))
601 icbm(ic_b, 3) = coset(lb(1), lb(2), lb(3) - 1)
602 END DO
603 DO ipgfb = 1, npgfb
604 zbx2 = zetb(ipgfb)*2.0_dp
605 icbx = (ipgfb - 1)*nb
606 DO ipgfa = 1, npgfa
607 icax = (ipgfa - 1)*nap
608 IF (.NOT. any(habd(icax + 1:icax + nap, icbx + 1:icbx + nb) /= 0.0_dp)) cycle
609 zax2 = zeta(ipgfa)*2.0_dp
610 icpa = (ipgfa - 1)*na
611 s0 = 0.0_dp
612 sp = 0.0_dp
613 sma = 0.0_dp
614 DO ic_b = ic_b_min, nb
615 lb(1:3) = indco(1:3, ic_b)
616 icob = icbx + ic_b
617 DO ic_a = ic_a_min, na
618 pval = pab(icpa + ic_a, icob)
619 s0 = s0 + pval*habd(icax + ic_a, icob)
620 sp(1) = sp(1) + pval*habd(icax + icap(ic_a, 1), icob)
621 sp(2) = sp(2) + pval*habd(icax + icap(ic_a, 2), icob)
622 sp(3) = sp(3) + pval*habd(icax + icap(ic_a, 3), icob)
623 sma(1) = sma(1) + pval*(rla(ic_a, 1)*habd(icax + icam(ic_a, 1), icob))
624 sma(2) = sma(2) + pval*(rla(ic_a, 2)*habd(icax + icam(ic_a, 2), icob))
625 sma(3) = sma(3) + pval*(rla(ic_a, 3)*habd(icax + icam(ic_a, 3), icob))
626 END DO
627 DO idir = 1, 3
628 IF (lb(idir) > 0) THEN
629 icob_m = icbx + icbm(ic_b, idir)
630 smb = 0.0_dp
631 DO ic_a = ic_a_min, na
632 smb = smb + pab(icpa + ic_a, icob)*habd(icax + ic_a, icob_m)
633 END DO
634 fb(idir) = fb(idir) - real(lb(idir), kind=dp)*smb
635 END IF
636 END DO
637 END DO
638 fa(1:3) = fa(1:3) + zax2*sp(1:3) - sma(1:3)
639 fb(1:3) = fb(1:3) + zbx2*(sp(1:3) - rab(1:3)*s0)
640 END DO
641 END DO
642
643 END SUBROUTINE verfc_force
644
645! **************************************************************************************************
646!> \brief Integrals = -Z*erfc(a*r)/r
647!> \param matrix_h ...
648!> \param matrix_p ...
649!> \param qs_kind_set ...
650!> \param atomic_kind_set ...
651!> \param particle_set ...
652!> \param calpha ...
653!> \param ccore ...
654!> \param eps_core_charge ...
655!> \param sab_orb ...
656!> \param sac_ae ...
657!> \param atcore ...
658! **************************************************************************************************
659 SUBROUTINE build_erfc(matrix_h, matrix_p, qs_kind_set, atomic_kind_set, particle_set, &
660 calpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore)
661
662 TYPE(dbcsr_p_type) :: matrix_h
663 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_p
664 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
665 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
666 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
667 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: calpha, ccore
668 REAL(kind=dp), INTENT(IN) :: eps_core_charge
669 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
670 POINTER :: sab_orb, sac_ae
671 REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
672 OPTIONAL :: atcore
673
674 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_erfc'
675
676 INTEGER :: handle, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, katom, kkind, &
677 ldai, ldsab, maxco, maxl, maxnset, maxsgf, mepos, na_plus, natom, nb_plus, ncoa, ncob, &
678 nij, nkind, nseta, nsetb, nthread, sgfa, sgfb, slot
679 INTEGER, DIMENSION(3) :: cellind
680 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
681 npgfb, nsgfa, nsgfb
682 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
683 LOGICAL :: doat, found
684 REAL(kind=dp) :: alpha_c, atk0, atk1, core_charge, &
685 core_radius, dab, dac, dbc, f0, rab2, &
686 rac2, rbc2, zeta_c
687 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ff
688 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
689 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
690 REAL(kind=dp), DIMENSION(3) :: rab, rac, rbc
691 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
692 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
693 sphi_b, zeta, zetb
694 TYPE(neighbor_list_iterator_p_type), &
695 DIMENSION(:), POINTER :: ap_iterator
696 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
697 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
698 TYPE(all_potential_type), POINTER :: all_potential
699 REAL(kind=dp), DIMENSION(SIZE(particle_set)) :: at_thread
700 TYPE(sgp_potential_type), POINTER :: sgp_potential
701
702!$ INTEGER(kind=omp_lock_kind), &
703!$ ALLOCATABLE, DIMENSION(:) :: locks
704!$ INTEGER :: lock_num, hash, hash1, hash2
705!$ INTEGER(KIND=int_8) :: iatom8
706!$ INTEGER, PARAMETER :: nlock = 501
707
708 mark_used(int_8)
709
710 CALL timeset(routinen, handle)
711
712 nkind = SIZE(atomic_kind_set)
713 natom = SIZE(particle_set)
714
715 doat = PRESENT(atcore)
716
717 IF (doat) THEN
718 IF (SIZE(matrix_p, 1) == 2) THEN
719 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
720 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
721 CALL dbcsr_add(matrix_p(2)%matrix, matrix_p(1)%matrix, &
722 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
723 END IF
724 END IF
725
726 at_thread = 0.0_dp
727
728 ALLOCATE (basis_set_list(nkind))
729 DO ikind = 1, nkind
730 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
731 IF (ASSOCIATED(basis_set_a)) THEN
732 basis_set_list(ikind)%gto_basis_set => basis_set_a
733 ELSE
734 NULLIFY (basis_set_list(ikind)%gto_basis_set)
735 END IF
736 END DO
737
738 CALL get_qs_kind_set(qs_kind_set, &
739 maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset)
740 CALL init_orbital_pointers(maxl + 1)
741 ldsab = max(maxco, maxsgf)
742 ldai = ncoset(maxl + 1)
743
744 nthread = 1
745!$ nthread = omp_get_max_threads()
746
747 ! iterator for basis/potential list
748 CALL neighbor_list_iterator_create(ap_iterator, sac_ae, search=.true., nthread=nthread)
749
750!$OMP PARALLEL &
751!$OMP DEFAULT (NONE) &
752!$OMP SHARED (ap_iterator, basis_set_list, &
753!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
754!$OMP sab_orb, sac_ae, nthread, ncoset, nkind, calpha, ccore, eps_core_charge, &
755!$OMP slot, ldsab, maxnset, ldai, maxl, maxco, doat, locks, natom) &
756!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
757!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, lock_num, &
758!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
759!$OMP zetb, zeta_c, alpha_c, core_charge, dab, irow, icol, h_block, found, iset, ncoa, &
760!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
761!$OMP rac, dac, rbc, rab2, rac2, rbc2, dbc, na_plus, nb_plus, verf, vnuc, &
762!$OMP set_radius_a, core_radius, rpgfa, mepos, &
763!$OMP atk0, atk1, habd, f0, katom, cellind, img, nij, ff, &
764!$OMP sgp_potential, all_potential, hash, hash1, hash2, iatom8) &
765!$OMP REDUCTION (+ : at_thread )
766
767!$OMP SINGLE
768!$ ALLOCATE (locks(nlock))
769!$OMP END SINGLE
770
771!$OMP DO
772!$ DO lock_num = 1, nlock
773!$ call omp_init_lock(locks(lock_num))
774!$ END DO
775!$OMP END DO
776
777 mepos = 0
778!$ mepos = omp_get_thread_num()
779
780 ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab))
781 ALLOCATE (verf(ldai, ldai, 2*maxl + 1), vnuc(ldai, ldai, 2*maxl + 1), ff(0:2*maxl))
782 IF (doat) THEN
783 ALLOCATE (pab(maxco, maxco, maxnset*maxnset))
784 END IF
785
786!$OMP DO SCHEDULE(GUIDED)
787 DO slot = 1, sab_orb(1)%nl_size
788
789 ikind = sab_orb(1)%nlist_task(slot)%ikind
790 jkind = sab_orb(1)%nlist_task(slot)%jkind
791 iatom = sab_orb(1)%nlist_task(slot)%iatom
792 jatom = sab_orb(1)%nlist_task(slot)%jatom
793 cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
794 rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
795
796 basis_set_a => basis_set_list(ikind)%gto_basis_set
797 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
798 basis_set_b => basis_set_list(jkind)%gto_basis_set
799 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
800!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
801!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
802 ! basis ikind
803 first_sgfa => basis_set_a%first_sgf
804 la_max => basis_set_a%lmax
805 la_min => basis_set_a%lmin
806 npgfa => basis_set_a%npgf
807 nseta = basis_set_a%nset
808 nsgfa => basis_set_a%nsgf_set
809 rpgfa => basis_set_a%pgf_radius
810 set_radius_a => basis_set_a%set_radius
811 sphi_a => basis_set_a%sphi
812 zeta => basis_set_a%zet
813 ! basis jkind
814 first_sgfb => basis_set_b%first_sgf
815 lb_max => basis_set_b%lmax
816 lb_min => basis_set_b%lmin
817 npgfb => basis_set_b%npgf
818 nsetb = basis_set_b%nset
819 nsgfb => basis_set_b%nsgf_set
820 rpgfb => basis_set_b%pgf_radius
821 set_radius_b => basis_set_b%set_radius
822 sphi_b => basis_set_b%sphi
823 zetb => basis_set_b%zet
824
825 dab = sqrt(sum(rab*rab))
826 img = 1
827
828 ! *** Use the symmetry of the first derivatives ***
829 IF (iatom == jatom) THEN
830 f0 = 1.0_dp
831 ELSE
832 f0 = 2.0_dp
833 END IF
834
835 ! *** Create matrix blocks for a new matrix block column ***
836 IF (iatom <= jatom) THEN
837 irow = iatom
838 icol = jatom
839 ELSE
840 irow = jatom
841 icol = iatom
842 END IF
843 NULLIFY (h_block)
844 CALL dbcsr_get_block_p(matrix=matrix_h%matrix, &
845 row=irow, col=icol, block=h_block, found=found)
846 IF (doat) THEN
847 NULLIFY (p_block)
848 CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
849 row=irow, col=icol, block=p_block, found=found)
850 cpassert(ASSOCIATED(p_block))
851 ! *** Decontract density matrix block ***
852 DO iset = 1, nseta
853 ncoa = npgfa(iset)*ncoset(la_max(iset))
854 sgfa = first_sgfa(1, iset)
855 DO jset = 1, nsetb
856 ncob = npgfb(jset)*ncoset(lb_max(jset))
857 sgfb = first_sgfb(1, jset)
858 nij = jset + (iset - 1)*maxnset
859 ! *** Decontract density matrix block ***
860 IF (iatom <= jatom) THEN
861 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
862 p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
863 ELSE
864 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
865 transpose(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1)))
866 END IF
867 pab(1:ncoa, 1:ncob, nij) = matmul(work(1:ncoa, 1:nsgfb(jset)), &
868 transpose(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
869 END DO
870 END DO
871 END IF
872
873 ! loop over all kinds of atoms
874 hab = 0._dp
875 DO kkind = 1, nkind
876 CALL get_qs_kind(qs_kind_set(kkind), zeff=zeta_c)
877 alpha_c = calpha(kkind)
878 core_charge = ccore(kkind)
879 core_radius = exp_radius(0, alpha_c, eps_core_charge, core_charge)
880 core_radius = max(core_radius, 10.0_dp)
881
882 CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
883
884 DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
885 CALL get_iterator_info(ap_iterator, jatom=katom, r=rac, mepos=mepos)
886
887 dac = sqrt(sum(rac*rac))
888 rbc(:) = rac(:) - rab(:)
889 dbc = sqrt(sum(rbc*rbc))
890 IF ((maxval(set_radius_a(:)) + core_radius < dac) .OR. &
891 (maxval(set_radius_b(:)) + core_radius < dbc)) THEN
892 cycle
893 END IF
894
895 DO iset = 1, nseta
896 IF (set_radius_a(iset) + core_radius < dac) cycle
897 ncoa = npgfa(iset)*ncoset(la_max(iset))
898 sgfa = first_sgfa(1, iset)
899 DO jset = 1, nsetb
900 IF (set_radius_b(jset) + core_radius < dbc) cycle
901 ncob = npgfb(jset)*ncoset(lb_max(jset))
902 sgfb = first_sgfb(1, jset)
903 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
904 rab2 = dab*dab
905 rac2 = dac*dac
906 rbc2 = dbc*dbc
907 nij = jset + (iset - 1)*maxnset
908 IF (doat) THEN
909 atk0 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
910 END IF
911 !
912 CALL verfc( &
913 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
914 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
915 alpha_c, core_radius, zeta_c, core_charge, &
916 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:))
917 !
918 IF (doat) THEN
919 atk1 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
920 at_thread(katom) = at_thread(katom) + (atk1 - atk0)
921 END IF
922 END DO
923 END DO
924 END DO
925 END DO
926 ! *** Contract nuclear attraction integrals
927 DO iset = 1, nseta
928 ncoa = npgfa(iset)*ncoset(la_max(iset))
929 sgfa = first_sgfa(1, iset)
930 DO jset = 1, nsetb
931 ncob = npgfb(jset)*ncoset(lb_max(jset))
932 sgfb = first_sgfb(1, jset)
933 nij = jset + (iset - 1)*maxnset
934!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
935!$ hash = MOD(hash1 + hash2, nlock) + 1
936 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab(1:ncoa, 1:ncob, nij), &
937 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
938!$ CALL omp_set_lock(locks(hash))
939 IF (iatom <= jatom) THEN
940 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
941 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
942 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
943 ELSE
944 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
945 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
946 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
947 END IF
948!$ CALL omp_unset_lock(locks(hash))
949 END DO
950 END DO
951
952 END DO
953
954 DEALLOCATE (hab, work, verf, vnuc, ff)
955
956!$OMP DO
957!$ DO lock_num = 1, nlock
958!$ call omp_destroy_lock(locks(lock_num))
959!$ END DO
960!$OMP END DO
961
962!$OMP SINGLE
963!$ DEALLOCATE (locks)
964!$OMP END SINGLE NOWAIT
965
966!$OMP END PARALLEL
967
968 CALL neighbor_list_iterator_release(ap_iterator)
969
970 DEALLOCATE (basis_set_list)
971
972 IF (doat) THEN
973 ! *** If LSD, then recover alpha density and beta density ***
974 ! *** from the total density (1) and the spin density (2) ***
975 IF (SIZE(matrix_p, 1) == 2) THEN
976 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
977 alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
978 CALL dbcsr_add(matrix_p(2)%matrix, matrix_p(1)%matrix, &
979 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
980 END IF
981 END IF
982
983 IF (doat) THEN
984 atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
985 END IF
986
987 CALL timestop(handle)
988
989 END SUBROUTINE build_erfc
990
991! **************************************************************************************************
992
993END MODULE core_ae
external_potential_types::get_potential
Definition external_potential_types.F:279
qs_neighbor_list_types::nl_sub_iterate
Definition qs_neighbor_list_types.F:96
ai_verfc
Build up the nuclear potential part of the core Hamiltonian matrix in the case of an allelectron calc...
Definition ai_verfc.F:66
ai_verfc::verfc
subroutine, public verfc(la_max1, npgfa, zeta, rpgfa, la_min1, lb_max1, npgfb, zetb, rpgfb, lb_min1, zetc, rpgfc, zc, cerf, rab, rab2, rac, rac2, rbc2, vabc, verf, vnuc, f, maxder, vabc_plus, vnabc, pvp_sum, pvp, dkh_erfc)
Calculation of the primitive three-center nuclear potential integrals <a|Z*erfc(r)/r|b> over Cartesia...
Definition ai_verfc.F:127
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
Definition ao_util.F:96
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:216
basis_set_types
Definition basis_set_types.F:15
core_ae
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
Definition core_ae.F:14
core_ae::build_core_ae
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index, basis_type, atcore)
...
Definition core_ae.F:87
core_ae::verfc_force
subroutine, public verfc_force(habd, pab, fa, fb, nder, la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, rab)
Forces from [a|V_erfc|b] nuclear attraction integrals for one shell pair.
Definition core_ae.F:560
core_ae::build_erfc
subroutine, public build_erfc(matrix_h, matrix_p, qs_kind_set, atomic_kind_set, particle_set, calpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore)
Integrals = -Z*erfc(a*r)/r.
Definition core_ae.F:661
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:483
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:976
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::nl_set_sub_iterator
subroutine, public nl_set_sub_iterator(iterator_set, ikind, jkind, iatom, mepos)
...
Definition qs_neighbor_list_types.F:288
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
external_potential_types::all_potential_type
Definition external_potential_types.F:75
external_potential_types::sgp_potential_type
Definition external_potential_types.F:179
particle_types::particle_type
Definition particle_types.F:35
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26