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