(git:ed6f26b)
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-2025 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
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 atcore ...
82! **************************************************************************************************
83 SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
84 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
85 nimages, cell_to_index, atcore)
86
87 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
88 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
89 TYPE(virial_type), POINTER :: virial
90 LOGICAL, INTENT(IN) :: calculate_forces
91 LOGICAL :: use_virial
92 INTEGER :: nder
93 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
94 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
95 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
96 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
97 POINTER :: sab_orb, sac_ae
98 INTEGER, INTENT(IN) :: nimages
99 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
100 REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
101 OPTIONAL :: atcore
102
103 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_core_ae'
104
105 INTEGER :: atom_a, handle, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, katom, &
106 kkind, ldai, ldsab, maxco, maxl, maxnset, maxsgf, mepos, na_plus, natom, nb_plus, ncoa, &
107 ncob, nij, nkind, nseta, nsetb, nthread, sgfa, sgfb, slot
108 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
109 INTEGER, DIMENSION(3) :: cellind
110 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
111 npgfb, nsgfa, nsgfb
112 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
113 LOGICAL :: doat, dokp, found
114 REAL(kind=dp) :: alpha_c, atk0, atk1, core_charge, &
115 core_radius, dab, dac, dbc, f0, rab2, &
116 rac2, rbc2, zeta_c
117 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ff
118 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
119 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
120 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
121 REAL(kind=dp), DIMENSION(3, 3) :: pv_thread
122 TYPE(neighbor_list_iterator_p_type), &
123 DIMENSION(:), POINTER :: ap_iterator
124 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
125 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
126 TYPE(all_potential_type), POINTER :: all_potential
127 REAL(kind=dp), DIMENSION(SIZE(particle_set)) :: at_thread
128 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
129 sphi_b, zeta, zetb
130 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
131 REAL(kind=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
132 TYPE(sgp_potential_type), POINTER :: sgp_potential
133
134!$ INTEGER(kind=omp_lock_kind), &
135!$ ALLOCATABLE, DIMENSION(:) :: locks
136!$ INTEGER :: lock_num, hash, hash1, hash2
137!$ INTEGER(KIND=int_8) :: iatom8
138!$ INTEGER, PARAMETER :: nlock = 501
139
140 mark_used(int_8)
141
142 IF (calculate_forces) THEN
143 CALL timeset(routinen//"_forces", handle)
144 ELSE
145 CALL timeset(routinen, handle)
146 END IF
147
148 nkind = SIZE(atomic_kind_set)
149 natom = SIZE(particle_set)
150
151 doat = PRESENT(atcore)
152 dokp = (nimages > 1)
153
154 IF (calculate_forces .OR. doat) THEN
155 IF (SIZE(matrix_p, 1) == 2) THEN
156 DO img = 1, nimages
157 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
158 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
159 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
160 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
161 END DO
162 END IF
163 END IF
164
165 force_thread = 0.0_dp
166 at_thread = 0.0_dp
167 pv_thread = 0.0_dp
168
169 ALLOCATE (basis_set_list(nkind))
170 DO ikind = 1, nkind
171 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
172 IF (ASSOCIATED(basis_set_a)) THEN
173 basis_set_list(ikind)%gto_basis_set => basis_set_a
174 ELSE
175 NULLIFY (basis_set_list(ikind)%gto_basis_set)
176 END IF
177 END DO
178
179 CALL get_qs_kind_set(qs_kind_set, &
180 maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset)
181 CALL init_orbital_pointers(maxl + nder + 1)
182 ldsab = max(maxco, maxsgf)
183 ldai = ncoset(maxl + nder + 1)
184
185 nthread = 1
186!$ nthread = omp_get_max_threads()
187
188 ! iterator for basis/potential list
189 CALL neighbor_list_iterator_create(ap_iterator, sac_ae, search=.true., nthread=nthread)
190
191!$OMP PARALLEL &
192!$OMP DEFAULT (NONE) &
193!$OMP SHARED (ap_iterator, basis_set_list, calculate_forces, use_virial, &
194!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
195!$OMP sab_orb, sac_ae, nthread, ncoset, nkind, cell_to_index, &
196!$OMP slot, ldsab, maxnset, ldai, nder, maxl, maxco, dokp, doat, locks, natom) &
197!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
198!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
199!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
200!$OMP zetb, zeta_c, alpha_c, core_charge, dab, irow, icol, h_block, found, iset, ncoa, &
201!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
202!$OMP rac, dac, rbc, rab2, rac2, rbc2, dbc, na_plus, nb_plus, verf, vnuc, &
203!$OMP set_radius_a, core_radius, rpgfa, force_a, force_b, mepos, &
204!$OMP atk0, atk1, habd, f0, katom, cellind, img, nij, ff, &
205!$OMP sgp_potential, all_potential, hash, hash1, hash2, iatom8) &
206!$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )
207
208!$OMP SINGLE
209!$ ALLOCATE (locks(nlock))
210!$OMP END SINGLE
211
212!$OMP DO
213!$ DO lock_num = 1, nlock
214!$ call omp_init_lock(locks(lock_num))
215!$ END DO
216!$OMP END DO
217
218 mepos = 0
219!$ mepos = omp_get_thread_num()
220
221 ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab))
222 ALLOCATE (verf(ldai, ldai, 2*maxl + nder + 1), vnuc(ldai, ldai, 2*maxl + nder + 1), ff(0:2*maxl + nder))
223 IF (calculate_forces .OR. doat) THEN
224 ALLOCATE (pab(maxco, maxco, maxnset*maxnset))
225 END IF
226
227!$OMP DO SCHEDULE(GUIDED)
228 DO slot = 1, sab_orb(1)%nl_size
229
230 ikind = sab_orb(1)%nlist_task(slot)%ikind
231 jkind = sab_orb(1)%nlist_task(slot)%jkind
232 iatom = sab_orb(1)%nlist_task(slot)%iatom
233 jatom = sab_orb(1)%nlist_task(slot)%jatom
234 cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
235 rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
236
237 basis_set_a => basis_set_list(ikind)%gto_basis_set
238 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
239 basis_set_b => basis_set_list(jkind)%gto_basis_set
240 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
241!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
242!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
243 ! basis ikind
244 first_sgfa => basis_set_a%first_sgf
245 la_max => basis_set_a%lmax
246 la_min => basis_set_a%lmin
247 npgfa => basis_set_a%npgf
248 nseta = basis_set_a%nset
249 nsgfa => basis_set_a%nsgf_set
250 rpgfa => basis_set_a%pgf_radius
251 set_radius_a => basis_set_a%set_radius
252 sphi_a => basis_set_a%sphi
253 zeta => basis_set_a%zet
254 ! basis jkind
255 first_sgfb => basis_set_b%first_sgf
256 lb_max => basis_set_b%lmax
257 lb_min => basis_set_b%lmin
258 npgfb => basis_set_b%npgf
259 nsetb = basis_set_b%nset
260 nsgfb => basis_set_b%nsgf_set
261 rpgfb => basis_set_b%pgf_radius
262 set_radius_b => basis_set_b%set_radius
263 sphi_b => basis_set_b%sphi
264 zetb => basis_set_b%zet
265
266 dab = sqrt(sum(rab*rab))
267
268 IF (dokp) THEN
269 img = cell_to_index(cellind(1), cellind(2), cellind(3))
270 ELSE
271 img = 1
272 END IF
273
274 ! *** Use the symmetry of the first derivatives ***
275 IF (iatom == jatom) THEN
276 f0 = 1.0_dp
277 ELSE
278 f0 = 2.0_dp
279 END IF
280
281 ! *** Create matrix blocks for a new matrix block column ***
282 IF (iatom <= jatom) THEN
283 irow = iatom
284 icol = jatom
285 ELSE
286 irow = jatom
287 icol = iatom
288 END IF
289 NULLIFY (h_block)
290 CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, &
291 row=irow, col=icol, block=h_block, found=found)
292 IF (calculate_forces .OR. doat) THEN
293 NULLIFY (p_block)
294 CALL dbcsr_get_block_p(matrix=matrix_p(1, img)%matrix, &
295 row=irow, col=icol, block=p_block, found=found)
296 cpassert(ASSOCIATED(p_block))
297 ! *** Decontract density matrix block ***
298 DO iset = 1, nseta
299 ncoa = npgfa(iset)*ncoset(la_max(iset))
300 sgfa = first_sgfa(1, iset)
301 DO jset = 1, nsetb
302 ncob = npgfb(jset)*ncoset(lb_max(jset))
303 sgfb = first_sgfb(1, jset)
304 nij = jset + (iset - 1)*maxnset
305 ! *** Decontract density matrix block ***
306 IF (iatom <= jatom) THEN
307 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
308 p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
309 ELSE
310 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
311 transpose(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1)))
312 END IF
313 pab(1:ncoa, 1:ncob, nij) = matmul(work(1:ncoa, 1:nsgfb(jset)), &
314 transpose(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
315 END DO
316 END DO
317 END IF
318
319 ! loop over all kinds for pseudopotential atoms
320 hab = 0._dp
321 DO kkind = 1, nkind
322 CALL get_qs_kind(qs_kind_set(kkind), all_potential=all_potential, &
323 sgp_potential=sgp_potential)
324 IF (ASSOCIATED(all_potential)) THEN
325 CALL get_potential(potential=all_potential, &
326 alpha_core_charge=alpha_c, zeff=zeta_c, &
327 ccore_charge=core_charge, core_charge_radius=core_radius)
328 ELSE IF (ASSOCIATED(sgp_potential)) THEN
329 CALL get_potential(potential=sgp_potential, &
330 alpha_core_charge=alpha_c, zeff=zeta_c, &
331 ccore_charge=core_charge, core_charge_radius=core_radius)
332 ELSE
333 cycle
334 END IF
335
336 CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
337
338 DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
339 CALL get_iterator_info(ap_iterator, jatom=katom, r=rac, mepos=mepos)
340
341 dac = sqrt(sum(rac*rac))
342 rbc(:) = rac(:) - rab(:)
343 dbc = sqrt(sum(rbc*rbc))
344 IF ((maxval(set_radius_a(:)) + core_radius < dac) .OR. &
345 (maxval(set_radius_b(:)) + core_radius < dbc)) THEN
346 cycle
347 END IF
348
349 DO iset = 1, nseta
350 IF (set_radius_a(iset) + core_radius < dac) cycle
351 ncoa = npgfa(iset)*ncoset(la_max(iset))
352 sgfa = first_sgfa(1, iset)
353 DO jset = 1, nsetb
354 IF (set_radius_b(jset) + core_radius < dbc) cycle
355 ncob = npgfb(jset)*ncoset(lb_max(jset))
356 sgfb = first_sgfb(1, jset)
357 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
358 rab2 = dab*dab
359 rac2 = dac*dac
360 rbc2 = dbc*dbc
361 nij = jset + (iset - 1)*maxnset
362 ! *** Calculate the GTH pseudo potential forces ***
363 IF (doat) THEN
364 atk0 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
365 END IF
366 IF (calculate_forces) THEN
367 na_plus = npgfa(iset)*ncoset(la_max(iset) + nder)
368 nb_plus = npgfb(jset)*ncoset(lb_max(jset))
369 ALLOCATE (habd(na_plus, nb_plus))
370 habd = 0._dp
371 CALL verfc( &
372 la_max(iset) + nder, npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
373 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
374 alpha_c, core_radius, zeta_c, core_charge, &
375 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:), &
376 nder, habd)
377
378 ! *** The derivatives w.r.t. atomic center c are ***
379 ! *** calculated using the translational invariance ***
380 ! *** of the first derivatives ***
381 CALL verfc_force(habd, pab(:, :, nij), force_a, force_b, nder, &
382 la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
383 lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), rab)
384
385 DEALLOCATE (habd)
386
387 force_thread(1, iatom) = force_thread(1, iatom) + f0*force_a(1)
388 force_thread(2, iatom) = force_thread(2, iatom) + f0*force_a(2)
389 force_thread(3, iatom) = force_thread(3, iatom) + f0*force_a(3)
390
391 force_thread(1, jatom) = force_thread(1, jatom) + f0*force_b(1)
392 force_thread(2, jatom) = force_thread(2, jatom) + f0*force_b(2)
393 force_thread(3, jatom) = force_thread(3, jatom) + f0*force_b(3)
394
395 force_thread(1, katom) = force_thread(1, katom) - f0*force_a(1) - f0*force_b(1)
396 force_thread(2, katom) = force_thread(2, katom) - f0*force_a(2) - f0*force_b(2)
397 force_thread(3, katom) = force_thread(3, katom) - f0*force_a(3) - f0*force_b(3)
398
399 IF (use_virial) THEN
400 CALL virial_pair_force(pv_thread, f0, force_a, rac)
401 CALL virial_pair_force(pv_thread, f0, force_b, rbc)
402 END IF
403 ELSE
404 CALL verfc( &
405 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
406 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
407 alpha_c, core_radius, zeta_c, core_charge, &
408 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:))
409 END IF
410 ! calculate atomic contributions
411 IF (doat) THEN
412 atk1 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
413 at_thread(katom) = at_thread(katom) + (atk1 - atk0)
414 END IF
415 END DO
416 END DO
417 END DO
418 END DO
419 ! *** Contract nuclear attraction integrals
420 DO iset = 1, nseta
421 ncoa = npgfa(iset)*ncoset(la_max(iset))
422 sgfa = first_sgfa(1, iset)
423 DO jset = 1, nsetb
424 ncob = npgfb(jset)*ncoset(lb_max(jset))
425 sgfb = first_sgfb(1, jset)
426 nij = jset + (iset - 1)*maxnset
427!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
428!$ hash = MOD(hash1 + hash2, nlock) + 1
429 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab(1:ncoa, 1:ncob, nij), &
430 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
431!$ CALL omp_set_lock(locks(hash))
432 IF (iatom <= jatom) THEN
433 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
434 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
435 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
436 ELSE
437 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
438 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
439 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
440 END IF
441!$ CALL omp_unset_lock(locks(hash))
442 END DO
443 END DO
444
445 END DO
446
447 DEALLOCATE (hab, work, verf, vnuc, ff)
448 IF (calculate_forces .OR. doat) THEN
449 DEALLOCATE (pab)
450 END IF
451
452!$OMP DO
453!$ DO lock_num = 1, nlock
454!$ call omp_destroy_lock(locks(lock_num))
455!$ END DO
456!$OMP END DO
457
458!$OMP SINGLE
459!$ DEALLOCATE (locks)
460!$OMP END SINGLE NOWAIT
461
462!$OMP END PARALLEL
463
464 CALL neighbor_list_iterator_release(ap_iterator)
465
466 DEALLOCATE (basis_set_list)
467
468 IF (calculate_forces .OR. doat) THEN
469 ! *** If LSD, then recover alpha density and beta density ***
470 ! *** from the total density (1) and the spin density (2) ***
471 IF (SIZE(matrix_p, 1) == 2) THEN
472 DO img = 1, nimages
473 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
474 alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
475 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
476 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
477 END DO
478 END IF
479 END IF
480
481 IF (calculate_forces) THEN
482 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
483!$OMP DO
484 DO iatom = 1, natom
485 atom_a = atom_of_kind(iatom)
486 ikind = kind_of(iatom)
487 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) + force_thread(:, iatom)
488 END DO
489!$OMP END DO
490 END IF
491 IF (doat) THEN
492 atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
493 END IF
494
495 IF (calculate_forces .AND. use_virial) THEN
496 virial%pv_ppl = virial%pv_ppl + pv_thread
497 virial%pv_virial = virial%pv_virial + pv_thread
498 END IF
499
500 CALL timestop(handle)
501
502 END SUBROUTINE build_core_ae
503
504! **************************************************************************************************
505!> \brief ...
506!> \param habd ...
507!> \param pab ...
508!> \param fa ...
509!> \param fb ...
510!> \param nder ...
511!> \param la_max ...
512!> \param la_min ...
513!> \param npgfa ...
514!> \param zeta ...
515!> \param lb_max ...
516!> \param lb_min ...
517!> \param npgfb ...
518!> \param zetb ...
519!> \param rab ...
520! **************************************************************************************************
521 SUBROUTINE verfc_force(habd, pab, fa, fb, nder, la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, rab)
522
523 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: habd, pab
524 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: fa, fb
525 INTEGER, INTENT(IN) :: nder, la_max, la_min, npgfa
526 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta
527 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
528 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zetb
529 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
530
531 INTEGER :: ic_a, ic_b, icam1, icam2, icam3, icap1, &
532 icap2, icap3, icax, icbm1, icbm2, &
533 icbm3, icbx, icoa, icob, ipgfa, ipgfb, &
534 na, nap, nb
535 INTEGER, DIMENSION(3) :: la, lb
536 REAL(kind=dp) :: zax2, zbx2
537
538 fa = 0.0_dp
539 fb = 0.0_dp
540
541 na = ncoset(la_max)
542 nap = ncoset(la_max + nder)
543 nb = ncoset(lb_max)
544 DO ipgfa = 1, npgfa
545 zax2 = zeta(ipgfa)*2.0_dp
546 DO ipgfb = 1, npgfb
547 zbx2 = zetb(ipgfb)*2.0_dp
548 DO ic_a = ncoset(la_min - 1) + 1, ncoset(la_max)
549 la(1:3) = indco(1:3, ic_a)
550 icap1 = coset(la(1) + 1, la(2), la(3))
551 icap2 = coset(la(1), la(2) + 1, la(3))
552 icap3 = coset(la(1), la(2), la(3) + 1)
553 icam1 = coset(la(1) - 1, la(2), la(3))
554 icam2 = coset(la(1), la(2) - 1, la(3))
555 icam3 = coset(la(1), la(2), la(3) - 1)
556 icoa = ic_a + (ipgfa - 1)*na
557 icax = (ipgfa - 1)*nap
558
559 DO ic_b = ncoset(lb_min - 1) + 1, ncoset(lb_max)
560 lb(1:3) = indco(1:3, ic_b)
561 icbm1 = coset(lb(1) - 1, lb(2), lb(3))
562 icbm2 = coset(lb(1), lb(2) - 1, lb(3))
563 icbm3 = coset(lb(1), lb(2), lb(3) - 1)
564 icob = ic_b + (ipgfb - 1)*nb
565 icbx = (ipgfb - 1)*nb
566
567 fa(1) = fa(1) - pab(icoa, icob)*(-zax2*habd(icap1 + icax, icob) + &
568 REAL(la(1), kind=dp)*habd(icam1 + icax, icob))
569 fa(2) = fa(2) - pab(icoa, icob)*(-zax2*habd(icap2 + icax, icob) + &
570 REAL(la(2), kind=dp)*habd(icam2 + icax, icob))
571 fa(3) = fa(3) - pab(icoa, icob)*(-zax2*habd(icap3 + icax, icob) + &
572 REAL(la(3), kind=dp)*habd(icam3 + icax, icob))
573
574 fb(1) = fb(1) - pab(icoa, icob)*( &
575 -zbx2*(habd(icap1 + icax, icob) - rab(1)*habd(ic_a + icax, icob)) + &
576 REAL(lb(1), kind=dp)*habd(ic_a + icax, icbm1 + icbx))
577 fb(2) = fb(2) - pab(icoa, icob)*( &
578 -zbx2*(habd(icap2 + icax, icob) - rab(2)*habd(ic_a + icax, icob)) + &
579 REAL(lb(2), kind=dp)*habd(ic_a + icax, icbm2 + icbx))
580 fb(3) = fb(3) - pab(icoa, icob)*( &
581 -zbx2*(habd(icap3 + icax, icob) - rab(3)*habd(ic_a + icax, icob)) + &
582 REAL(lb(3), kind=dp)*habd(ic_a + icax, icbm3 + icbx))
583
584 END DO ! ic_b
585 END DO ! ic_a
586 END DO ! ipgfb
587 END DO ! ipgfa
588
589 END SUBROUTINE verfc_force
590
591! **************************************************************************************************
592!> \brief Integrals = -Z*erfc(a*r)/r
593!> \param matrix_h ...
594!> \param matrix_p ...
595!> \param qs_kind_set ...
596!> \param atomic_kind_set ...
597!> \param particle_set ...
598!> \param calpha ...
599!> \param ccore ...
600!> \param eps_core_charge ...
601!> \param sab_orb ...
602!> \param sac_ae ...
603!> \param atcore ...
604! **************************************************************************************************
605 SUBROUTINE build_erfc(matrix_h, matrix_p, qs_kind_set, atomic_kind_set, particle_set, &
606 calpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore)
607
608 TYPE(dbcsr_p_type) :: matrix_h
609 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_p
610 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
611 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
612 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
613 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: calpha, ccore
614 REAL(kind=dp), INTENT(IN) :: eps_core_charge
615 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
616 POINTER :: sab_orb, sac_ae
617 REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
618 OPTIONAL :: atcore
619
620 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_erfc'
621
622 INTEGER :: handle, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, katom, kkind, &
623 ldai, ldsab, maxco, maxl, maxnset, maxsgf, mepos, na_plus, natom, nb_plus, ncoa, ncob, &
624 nij, nkind, nseta, nsetb, nthread, sgfa, sgfb, slot
625 INTEGER, DIMENSION(3) :: cellind
626 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
627 npgfb, nsgfa, nsgfb
628 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
629 LOGICAL :: doat, found
630 REAL(kind=dp) :: alpha_c, atk0, atk1, core_charge, &
631 core_radius, dab, dac, dbc, f0, rab2, &
632 rac2, rbc2, zeta_c
633 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ff
634 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
635 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
636 REAL(kind=dp), DIMENSION(3) :: rab, rac, rbc
637 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
638 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
639 sphi_b, zeta, zetb
640 TYPE(neighbor_list_iterator_p_type), &
641 DIMENSION(:), POINTER :: ap_iterator
642 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
643 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
644 TYPE(all_potential_type), POINTER :: all_potential
645 REAL(kind=dp), DIMENSION(SIZE(particle_set)) :: at_thread
646 TYPE(sgp_potential_type), POINTER :: sgp_potential
647
648!$ INTEGER(kind=omp_lock_kind), &
649!$ ALLOCATABLE, DIMENSION(:) :: locks
650!$ INTEGER :: lock_num, hash, hash1, hash2
651!$ INTEGER(KIND=int_8) :: iatom8
652!$ INTEGER, PARAMETER :: nlock = 501
653
654 mark_used(int_8)
655
656 CALL timeset(routinen, handle)
657
658 nkind = SIZE(atomic_kind_set)
659 natom = SIZE(particle_set)
660
661 doat = PRESENT(atcore)
662
663 IF (doat) THEN
664 IF (SIZE(matrix_p, 1) == 2) THEN
665 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
666 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
667 CALL dbcsr_add(matrix_p(2)%matrix, matrix_p(1)%matrix, &
668 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
669 END IF
670 END IF
671
672 at_thread = 0.0_dp
673
674 ALLOCATE (basis_set_list(nkind))
675 DO ikind = 1, nkind
676 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
677 IF (ASSOCIATED(basis_set_a)) THEN
678 basis_set_list(ikind)%gto_basis_set => basis_set_a
679 ELSE
680 NULLIFY (basis_set_list(ikind)%gto_basis_set)
681 END IF
682 END DO
683
684 CALL get_qs_kind_set(qs_kind_set, &
685 maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset)
686 CALL init_orbital_pointers(maxl + 1)
687 ldsab = max(maxco, maxsgf)
688 ldai = ncoset(maxl + 1)
689
690 nthread = 1
691!$ nthread = omp_get_max_threads()
692
693 ! iterator for basis/potential list
694 CALL neighbor_list_iterator_create(ap_iterator, sac_ae, search=.true., nthread=nthread)
695
696!$OMP PARALLEL &
697!$OMP DEFAULT (NONE) &
698!$OMP SHARED (ap_iterator, basis_set_list, &
699!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
700!$OMP sab_orb, sac_ae, nthread, ncoset, nkind, calpha, ccore, eps_core_charge, &
701!$OMP slot, ldsab, maxnset, ldai, maxl, maxco, doat, locks, natom) &
702!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
703!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
704!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
705!$OMP zetb, zeta_c, alpha_c, core_charge, dab, irow, icol, h_block, found, iset, ncoa, &
706!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
707!$OMP rac, dac, rbc, rab2, rac2, rbc2, dbc, na_plus, nb_plus, verf, vnuc, &
708!$OMP set_radius_a, core_radius, rpgfa, mepos, &
709!$OMP atk0, atk1, habd, f0, katom, cellind, img, nij, ff, &
710!$OMP sgp_potential, all_potential, hash, hash1, hash2, iatom8) &
711!$OMP REDUCTION (+ : at_thread )
712
713!$OMP SINGLE
714!$ ALLOCATE (locks(nlock))
715!$OMP END SINGLE
716
717!$OMP DO
718!$ DO lock_num = 1, nlock
719!$ call omp_init_lock(locks(lock_num))
720!$ END DO
721!$OMP END DO
722
723 mepos = 0
724!$ mepos = omp_get_thread_num()
725
726 ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab))
727 ALLOCATE (verf(ldai, ldai, 2*maxl + 1), vnuc(ldai, ldai, 2*maxl + 1), ff(0:2*maxl))
728 IF (doat) THEN
729 ALLOCATE (pab(maxco, maxco, maxnset*maxnset))
730 END IF
731
732!$OMP DO SCHEDULE(GUIDED)
733 DO slot = 1, sab_orb(1)%nl_size
734
735 ikind = sab_orb(1)%nlist_task(slot)%ikind
736 jkind = sab_orb(1)%nlist_task(slot)%jkind
737 iatom = sab_orb(1)%nlist_task(slot)%iatom
738 jatom = sab_orb(1)%nlist_task(slot)%jatom
739 cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
740 rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
741
742 basis_set_a => basis_set_list(ikind)%gto_basis_set
743 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
744 basis_set_b => basis_set_list(jkind)%gto_basis_set
745 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
746!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
747!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
748 ! basis ikind
749 first_sgfa => basis_set_a%first_sgf
750 la_max => basis_set_a%lmax
751 la_min => basis_set_a%lmin
752 npgfa => basis_set_a%npgf
753 nseta = basis_set_a%nset
754 nsgfa => basis_set_a%nsgf_set
755 rpgfa => basis_set_a%pgf_radius
756 set_radius_a => basis_set_a%set_radius
757 sphi_a => basis_set_a%sphi
758 zeta => basis_set_a%zet
759 ! basis jkind
760 first_sgfb => basis_set_b%first_sgf
761 lb_max => basis_set_b%lmax
762 lb_min => basis_set_b%lmin
763 npgfb => basis_set_b%npgf
764 nsetb = basis_set_b%nset
765 nsgfb => basis_set_b%nsgf_set
766 rpgfb => basis_set_b%pgf_radius
767 set_radius_b => basis_set_b%set_radius
768 sphi_b => basis_set_b%sphi
769 zetb => basis_set_b%zet
770
771 dab = sqrt(sum(rab*rab))
772 img = 1
773
774 ! *** Use the symmetry of the first derivatives ***
775 IF (iatom == jatom) THEN
776 f0 = 1.0_dp
777 ELSE
778 f0 = 2.0_dp
779 END IF
780
781 ! *** Create matrix blocks for a new matrix block column ***
782 IF (iatom <= jatom) THEN
783 irow = iatom
784 icol = jatom
785 ELSE
786 irow = jatom
787 icol = iatom
788 END IF
789 NULLIFY (h_block)
790 CALL dbcsr_get_block_p(matrix=matrix_h%matrix, &
791 row=irow, col=icol, block=h_block, found=found)
792 IF (doat) THEN
793 NULLIFY (p_block)
794 CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
795 row=irow, col=icol, block=p_block, found=found)
796 cpassert(ASSOCIATED(p_block))
797 ! *** Decontract density matrix block ***
798 DO iset = 1, nseta
799 ncoa = npgfa(iset)*ncoset(la_max(iset))
800 sgfa = first_sgfa(1, iset)
801 DO jset = 1, nsetb
802 ncob = npgfb(jset)*ncoset(lb_max(jset))
803 sgfb = first_sgfb(1, jset)
804 nij = jset + (iset - 1)*maxnset
805 ! *** Decontract density matrix block ***
806 IF (iatom <= jatom) THEN
807 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
808 p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
809 ELSE
810 work(1:ncoa, 1:nsgfb(jset)) = matmul(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
811 transpose(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1)))
812 END IF
813 pab(1:ncoa, 1:ncob, nij) = matmul(work(1:ncoa, 1:nsgfb(jset)), &
814 transpose(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
815 END DO
816 END DO
817 END IF
818
819 ! loop over all kinds of atoms
820 hab = 0._dp
821 DO kkind = 1, nkind
822 CALL get_qs_kind(qs_kind_set(kkind), zeff=zeta_c)
823 alpha_c = calpha(kkind)
824 core_charge = ccore(kkind)
825 core_radius = exp_radius(0, alpha_c, eps_core_charge, core_charge)
826 core_radius = max(core_radius, 10.0_dp)
827
828 CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
829
830 DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
831 CALL get_iterator_info(ap_iterator, jatom=katom, r=rac, mepos=mepos)
832
833 dac = sqrt(sum(rac*rac))
834 rbc(:) = rac(:) - rab(:)
835 dbc = sqrt(sum(rbc*rbc))
836 IF ((maxval(set_radius_a(:)) + core_radius < dac) .OR. &
837 (maxval(set_radius_b(:)) + core_radius < dbc)) THEN
838 cycle
839 END IF
840
841 DO iset = 1, nseta
842 IF (set_radius_a(iset) + core_radius < dac) cycle
843 ncoa = npgfa(iset)*ncoset(la_max(iset))
844 sgfa = first_sgfa(1, iset)
845 DO jset = 1, nsetb
846 IF (set_radius_b(jset) + core_radius < dbc) cycle
847 ncob = npgfb(jset)*ncoset(lb_max(jset))
848 sgfb = first_sgfb(1, jset)
849 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
850 rab2 = dab*dab
851 rac2 = dac*dac
852 rbc2 = dbc*dbc
853 nij = jset + (iset - 1)*maxnset
854 IF (doat) THEN
855 atk0 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
856 END IF
857 !
858 CALL verfc( &
859 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
860 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
861 alpha_c, core_radius, zeta_c, core_charge, &
862 rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:))
863 !
864 IF (doat) THEN
865 atk1 = f0*sum(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
866 at_thread(katom) = at_thread(katom) + (atk1 - atk0)
867 END IF
868 END DO
869 END DO
870 END DO
871 END DO
872 ! *** Contract nuclear attraction integrals
873 DO iset = 1, nseta
874 ncoa = npgfa(iset)*ncoset(la_max(iset))
875 sgfa = first_sgfa(1, iset)
876 DO jset = 1, nsetb
877 ncob = npgfb(jset)*ncoset(lb_max(jset))
878 sgfb = first_sgfb(1, jset)
879 nij = jset + (iset - 1)*maxnset
880!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
881!$ hash = MOD(hash1 + hash2, nlock) + 1
882 work(1:ncoa, 1:nsgfb(jset)) = matmul(hab(1:ncoa, 1:ncob, nij), &
883 sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
884!$ CALL omp_set_lock(locks(hash))
885 IF (iatom <= jatom) THEN
886 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
887 h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
888 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
889 ELSE
890 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
891 h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
892 matmul(transpose(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
893 END IF
894!$ CALL omp_unset_lock(locks(hash))
895 END DO
896 END DO
897
898 END DO
899
900 DEALLOCATE (hab, work, verf, vnuc, ff)
901
902!$OMP DO
903!$ DO lock_num = 1, nlock
904!$ call omp_destroy_lock(locks(lock_num))
905!$ END DO
906!$OMP END DO
907
908!$OMP SINGLE
909!$ DEALLOCATE (locks)
910!$OMP END SINGLE NOWAIT
911
912!$OMP END PARALLEL
913
914 CALL neighbor_list_iterator_release(ap_iterator)
915
916 DEALLOCATE (basis_set_list)
917
918 IF (doat) THEN
919 ! *** If LSD, then recover alpha density and beta density ***
920 ! *** from the total density (1) and the spin density (2) ***
921 IF (SIZE(matrix_p, 1) == 2) THEN
922 CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
923 alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
924 CALL dbcsr_add(matrix_p(2)%matrix, matrix_p(1)%matrix, &
925 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
926 END IF
927 END IF
928
929 IF (doat) THEN
930 atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
931 END IF
932
933 CALL timestop(handle)
934
935 END SUBROUTINE build_erfc
936
937! **************************************************************************************************
938
939END MODULE core_ae
external_potential_types::get_potential
Definition external_potential_types.F:276
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:82
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:141
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:223
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, atcore)
...
Definition core_ae.F:86
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:607
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:679
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
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
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, 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, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
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:549
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F: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:170
external_potential_types::all_potential_type
Definition external_potential_types.F:70
external_potential_types::sgp_potential_type
Definition external_potential_types.F:172
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:171
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