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