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