(git:1f285aa)
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 ! **************************************************************************************************
14 MODULE core_ae
15  USE ai_verfc, ONLY: verfc
16  USE ao_util, ONLY: exp_radius
17  USE atomic_kind_types, ONLY: atomic_kind_type,&
19  USE basis_set_types, ONLY: gto_basis_set_p_type,&
20  gto_basis_set_type
21  USE 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,&
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,&
37  qs_kind_type
40  neighbor_list_iterator_p_type,&
42  neighbor_list_set_p_type,&
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 
63 CONTAINS
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
117  TYPE(neighbor_list_iterator_p_type), &
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
618  TYPE(neighbor_list_iterator_p_type), &
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 
847 END 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.
Definition: qs_kind_types.F:23
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.
integer function, public nl_sub_iterate(iterator_set, mepos)
...
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.