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