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