Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : MODULE qs_vcd_ao
8 :
9 : USE ai_cartesian_operators, ONLY: apply_cartesian_operator
10 : USE ai_contraction, ONLY: block_add,&
11 : contraction
12 : USE ai_kinetic, ONLY: kinetic
13 : USE ai_overlap_ppl, ONLY: ppl_integral
14 : USE ao_util, ONLY: exp_radius_very_extended
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind_set
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_p_type,&
19 : gto_basis_set_type
20 : USE block_p_types, ONLY: block_p_type
21 : USE cell_types, ONLY: cell_type,&
22 : pbc
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_add, dbcsr_copy, dbcsr_desymmetrize, dbcsr_distribution_get, &
26 : dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_init_p, &
27 : dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_work_create
28 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
29 : dbcsr_deallocate_matrix_set
30 : USE external_potential_types, ONLY: get_potential,&
31 : gth_potential_type,&
32 : sgp_potential_type
33 : USE gaussian_gridlevels, ONLY: gridlevel_info_type
34 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 : section_vals_type
36 : USE kinds, ONLY: default_string_length,&
37 : dp,&
38 : int_8
39 : USE memory_utilities, ONLY: reallocate
40 : USE message_passing, ONLY: mp_comm_type
41 : USE orbital_pointers, ONLY: init_orbital_pointers,&
42 : ncoset
43 : USE particle_types, ONLY: particle_type
44 : USE pw_env_types, ONLY: pw_env_get,&
45 : pw_env_type
46 : USE pw_methods, ONLY: pw_axpy,&
47 : pw_zero
48 : USE pw_pool_types, ONLY: pw_pool_type
49 : USE pw_types, ONLY: pw_r3d_rs_type
50 : USE qs_energy_types, ONLY: qs_energy_type
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
54 : get_memory_usage
55 : USE qs_integrate_potential, ONLY: integrate_pgf_product
56 : USE qs_kind_types, ONLY: get_qs_kind,&
57 : get_qs_kind_set,&
58 : qs_kind_type
59 : USE qs_ks_types, ONLY: qs_ks_env_type
60 : USE qs_linres_types, ONLY: vcd_env_type
61 : USE qs_moments, ONLY: build_local_moment_matrix
62 : USE qs_neighbor_list_types, ONLY: &
63 : get_iterator_info, get_neighbor_list_set_p, neighbor_list_iterate, &
64 : neighbor_list_iterator_create, neighbor_list_iterator_p_type, &
65 : neighbor_list_iterator_release, neighbor_list_set_p_type, nl_set_sub_iterator, &
66 : nl_sub_iterate
67 : USE qs_rho_types, ONLY: qs_rho_type
68 : USE qs_vxc, ONLY: qs_vxc_create
69 : USE realspace_grid_types, ONLY: realspace_grid_type
70 : USE rs_pw_interface, ONLY: potential_pw2rs
71 : USE sap_kind_types, ONLY: alist_type,&
72 : build_sap_ints,&
73 : get_alist,&
74 : release_sap_int,&
75 : sap_int_type,&
76 : sap_sort
77 : USE task_list_types, ONLY: atom_pair_type,&
78 : task_list_type,&
79 : task_type
80 : USE upf_local_integrals, ONLY: upf_local_init,&
81 : upf_local_matrix
82 : USE virial_types, ONLY: virial_type
83 :
84 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
85 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
86 : !$ omp_init_lock, omp_set_lock, &
87 : !$ omp_unset_lock, omp_destroy_lock
88 :
89 : #include "./base/base_uses.f90"
90 :
91 : IMPLICIT NONE
92 :
93 : PRIVATE
94 :
95 : ! *** Global parameters ***
96 :
97 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd_ao'
98 : INTEGER, PARAMETER :: bi_1 = 1, bi_x = 2, bi_y = 3, bi_z = 4, bi_xx = 5, &
99 : bi_xy = 6, bi_xz = 7, bi_yy = 8, bi_yz = 9, bi_zz = 10
100 : INTEGER, DIMENSION(3), PARAMETER :: bi_r = [bi_x, bi_y, bi_z]
101 : INTEGER, DIMENSION(3, 3), PARAMETER :: bi_rr = RESHAPE([bi_xx, bi_xy, bi_xz, bi_xy, bi_yy, bi_yz, &
102 : bi_xz, bi_yz, bi_zz], [3, 3])
103 :
104 : ! *** Public subroutines ***
105 :
106 : PUBLIC :: build_dSdV_matrix, build_com_rpnl_r, &
107 : hr_mult_by_delta_3d, build_dcom_rpnl, &
108 : build_matrix_hr_rh
109 :
110 : CONTAINS
111 :
112 : ! **************************************************************************************************
113 : !> \brief Build the matrix Hr*delta_nu^\lambda - rH*delta_mu^\lambda
114 : !> \param vcd_env ...
115 : !> \param qs_env ...
116 : !> \param rc ...
117 : !> \author Edward Ditler
118 : ! **************************************************************************************************
119 38 : SUBROUTINE build_matrix_hr_rh(vcd_env, qs_env, rc)
120 : TYPE(vcd_env_type) :: vcd_env
121 : TYPE(qs_environment_type), POINTER :: qs_env
122 : REAL(dp), DIMENSION(3) :: rc
123 :
124 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_matrix_hr_rh'
125 : INTEGER, PARAMETER :: ispin = 1
126 :
127 : INTEGER :: handle, i
128 : TYPE(cell_type), POINTER :: cell
129 : TYPE(dft_control_type), POINTER :: dft_control
130 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
131 38 : POINTER :: sab_all, sap_ppnl
132 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
133 38 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
134 :
135 38 : CALL timeset(routineN, handle)
136 :
137 : CALL get_qs_env(qs_env=qs_env, &
138 : dft_control=dft_control, &
139 : particle_set=particle_set, &
140 : sab_all=sab_all, &
141 : sap_ppnl=sap_ppnl, &
142 : qs_kind_set=qs_kind_set, &
143 38 : cell=cell)
144 :
145 152 : DO i = 1, 3
146 114 : CALL dbcsr_set(vcd_env%matrix_hr(ispin, i)%matrix, 0._dp)
147 152 : CALL dbcsr_set(vcd_env%matrix_rh(ispin, i)%matrix, 0._dp)
148 : END DO
149 :
150 : ASSOCIATE (matrix_hr_1d => vcd_env%matrix_hr(ispin, 1:3), &
151 : matrix_rh_1d => vcd_env%matrix_rh(ispin, 1:3))
152 : CALL build_rpnl_matrices(matrix_hr_1d, matrix_rh_1d, qs_kind_set, particle_set, sab_all, sap_ppnl, &
153 38 : dft_control%qs_control%eps_ppnl, cell, rc)
154 38 : CALL build_tr_matrices(matrix_hr_1d, matrix_rh_1d, qs_env, qs_kind_set, "ORB", sab_all, rc)
155 38 : CALL build_rcore_matrices(matrix_hr_1d, matrix_rh_1d, qs_env, qs_kind_set, "ORB", sab_all, rc)
156 76 : CALL build_vhxc_matrices(vcd_env%matrix_hr, vcd_env%matrix_rh, qs_env, rc)
157 : END ASSOCIATE
158 :
159 38 : CALL timestop(handle)
160 38 : END SUBROUTINE build_matrix_hr_rh
161 :
162 : ! **************************************************************************************************
163 : !> \brief Build the ordered products V_nl*r and r*V_nl.
164 : !> \param matrix_vr ...
165 : !> \param matrix_rv ...
166 : !> \param qs_kind_set ...
167 : !> \param particle_set ...
168 : !> \param sab_all ...
169 : !> \param sap_ppnl ...
170 : !> \param eps_ppnl ...
171 : !> \param cell ...
172 : !> \param ref_point ...
173 : !> \author Edward Ditler, Tomas Zimmermann
174 : ! **************************************************************************************************
175 38 : SUBROUTINE build_rpnl_matrices(matrix_vr, matrix_rv, qs_kind_set, particle_set, sab_all, sap_ppnl, eps_ppnl, &
176 : cell, ref_point)
177 :
178 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_vr, matrix_rv
179 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
180 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
181 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
182 : POINTER :: sab_all, sap_ppnl
183 : REAL(KIND=dp), INTENT(IN) :: eps_ppnl
184 : TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
185 : REAL(KIND=dp), DIMENSION(3) :: ref_point
186 :
187 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_rpnl_matrices'
188 :
189 : INTEGER :: handle, i, iab, iac, iatom, ibc, icol, &
190 : ikind, irow, jatom, jkind, kac, kbc, &
191 : kkind, na, natom, nb, nkind, np, slot
192 : INTEGER, DIMENSION(3) :: cell_b
193 : LOGICAL :: found, ppnl_present
194 38 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
195 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
196 266 : TYPE(block_p_type), DIMENSION(3) :: blocks_rv, blocks_vr
197 38 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
198 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
199 38 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
200 :
201 : !$ INTEGER(kind=omp_lock_kind), &
202 38 : !$ ALLOCATABLE, DIMENSION(:) :: locks
203 : !$ INTEGER :: lock_num, hash
204 : !$ INTEGER, PARAMETER :: nlock = 501
205 :
206 38 : ppnl_present = ASSOCIATED(sap_ppnl)
207 38 : IF (.NOT. ppnl_present) RETURN
208 :
209 38 : CALL timeset(routineN, handle)
210 38 : nkind = SIZE(qs_kind_set)
211 38 : natom = SIZE(particle_set)
212 :
213 : ! sap_int needs to be shared as multiple threads need to access this
214 38 : NULLIFY (sap_int)
215 266 : ALLOCATE (sap_int(nkind*nkind))
216 190 : DO i = 1, nkind*nkind
217 152 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
218 190 : sap_int(i)%nalist = 0
219 : END DO
220 :
221 : MARK_USED(ref_point)
222 : ! "nder" in moment_mode is "order"
223 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder=1, moment_mode=.TRUE., &
224 38 : particle_set=particle_set, cell=cell, refpoint=ref_point)
225 :
226 : ! *** Set up a sorting index
227 38 : CALL sap_sort(sap_int)
228 :
229 190 : ALLOCATE (basis_set(nkind))
230 114 : DO ikind = 1, nkind
231 76 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
232 114 : IF (ASSOCIATED(orb_basis_set)) THEN
233 76 : basis_set(ikind)%gto_basis_set => orb_basis_set
234 : ELSE
235 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
236 : END IF
237 : END DO
238 :
239 : ! *** All integrals needed have been calculated and stored in sap_int
240 : ! *** We now calculate the ordered product matrix elements
241 :
242 : !$OMP PARALLEL &
243 : !$OMP DEFAULT (NONE) &
244 : !$OMP SHARED (basis_set, matrix_vr, matrix_rv, &
245 : !$OMP sap_int, nkind, eps_ppnl, locks, sab_all) &
246 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, &
247 : !$OMP iab, irow, icol, blocks_vr, blocks_rv, &
248 : !$OMP found, iac, ibc, alist_ac, alist_bc, &
249 : !$OMP na, np, nb, kkind, kac, kbc, i, lock_num, &
250 38 : !$OMP hash, natom, acint, bcint, achint, bchint)
251 :
252 : !$OMP SINGLE
253 : !$ ALLOCATE (locks(nlock))
254 : !$OMP END SINGLE
255 :
256 : !$OMP DO
257 : !$ DO lock_num = 1, nlock
258 : !$ call omp_init_lock(locks(lock_num))
259 : !$ END DO
260 : !$OMP END DO
261 :
262 : !$OMP DO SCHEDULE(GUIDED)
263 :
264 : DO slot = 1, sab_all(1)%nl_size
265 :
266 : ikind = sab_all(1)%nlist_task(slot)%ikind
267 : jkind = sab_all(1)%nlist_task(slot)%jkind
268 : iatom = sab_all(1)%nlist_task(slot)%iatom
269 : jatom = sab_all(1)%nlist_task(slot)%jatom
270 : cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
271 :
272 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
273 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
274 : iab = ikind + nkind*(jkind - 1)
275 :
276 : ! *** Create matrix blocks for a new matrix block column ***
277 : irow = iatom
278 : icol = jatom
279 : DO i = 1, 3
280 : CALL dbcsr_get_block_p(matrix_vr(i)%matrix, irow, icol, blocks_vr(i)%block, found)
281 : CPASSERT(found)
282 : CALL dbcsr_get_block_p(matrix_rv(i)%matrix, irow, icol, blocks_rv(i)%block, found)
283 : CPASSERT(found)
284 : END DO
285 :
286 : ! loop over all kinds for projector atom
287 : DO kkind = 1, nkind
288 : iac = ikind + nkind*(kkind - 1)
289 : ibc = jkind + nkind*(kkind - 1)
290 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
291 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
292 : CALL get_alist(sap_int(iac), alist_ac, iatom)
293 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
294 :
295 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
296 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
297 : DO kac = 1, alist_ac%nclist
298 : DO kbc = 1, alist_bc%nclist
299 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
300 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
301 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
302 : acint => alist_ac%clist(kac)%acint
303 : bcint => alist_bc%clist(kbc)%acint
304 : achint => alist_ac%clist(kac)%achint
305 : bchint => alist_bc%clist(kbc)%achint
306 : na = SIZE(acint, 1)
307 : np = SIZE(acint, 2)
308 : nb = SIZE(bcint, 1)
309 : !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
310 : !$ CALL omp_set_lock(locks(hash))
311 : DO i = 1, 3
312 : ! Vnl*r
313 : blocks_vr(i)%block(1:na, 1:nb) = blocks_vr(i)%block(1:na, 1:nb) + &
314 : MATMUL(achint(1:na, 1:np, 1), &
315 : TRANSPOSE(bcint(1:nb, 1:np, i + 1)))
316 : ! r*Vnl
317 : blocks_rv(i)%block(1:na, 1:nb) = blocks_rv(i)%block(1:na, 1:nb) + &
318 : MATMUL(achint(1:na, 1:np, i + 1), &
319 : TRANSPOSE(bcint(1:nb, 1:np, 1)))
320 : END DO
321 : !$ CALL omp_unset_lock(locks(hash))
322 : EXIT ! We have found a match and there can be only one single match
323 : END IF
324 : END DO
325 : END DO
326 : END DO
327 : DO i = 1, 3
328 : NULLIFY (blocks_vr(i)%block)
329 : NULLIFY (blocks_rv(i)%block)
330 : END DO
331 : END DO
332 :
333 : !$OMP DO
334 : !$ DO lock_num = 1, nlock
335 : !$ call omp_destroy_lock(locks(lock_num))
336 : !$ END DO
337 : !$OMP END DO
338 :
339 : !$OMP SINGLE
340 : !$ DEALLOCATE (locks)
341 : !$OMP END SINGLE NOWAIT
342 :
343 : !$OMP END PARALLEL
344 :
345 38 : CALL release_sap_int(sap_int)
346 :
347 38 : DEALLOCATE (basis_set)
348 :
349 38 : CALL timestop(handle)
350 :
351 114 : END SUBROUTINE build_rpnl_matrices
352 :
353 : ! **************************************************************************************************
354 : !> \brief Calculation of the product Tr or rT over Cartesian Gaussian functions.
355 : !> \param matrix_tr ...
356 : !> \param matrix_rt ...
357 : !> \param qs_env ...
358 : !> \param qs_kind_set ...
359 : !> \param basis_type basis set to be used
360 : !> \param sab_nl pair list (must be consistent with basis sets!)
361 : !> \param rc ...
362 : !> \date 11.10.2010
363 : !> \par History
364 : !> Ported from qs_overlap, replaces code in build_core_hamiltonian
365 : !> Refactoring [07.2014] JGH
366 : !> Simplify options and use new kinetic energy integral routine
367 : !> Adapted from qs_kinetic [07.2016]
368 : !> Adapted from the kinetic position-product recurrence [2021] by ED
369 : !> \author JGH
370 : !> \version 1.0
371 : ! **************************************************************************************************
372 38 : SUBROUTINE build_tr_matrices(matrix_tr, matrix_rt, qs_env, qs_kind_set, basis_type, sab_nl, rc)
373 :
374 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_tr, matrix_rt
375 : TYPE(qs_environment_type), POINTER :: qs_env
376 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
377 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
378 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
379 : POINTER :: sab_nl
380 : REAL(KIND=dp), DIMENSION(3) :: rc
381 :
382 : CHARACTER(len=*), PARAMETER :: routineN = 'build_tr_matrices'
383 :
384 : INTEGER :: handle, i, iatom, icol, ikind, ir, irow, &
385 : iset, jatom, jkind, jset, ldsab, ltab, &
386 : natom, ncoa, ncob, nkind, nseta, &
387 : nsetb, sgfa, sgfb, slot
388 : INTEGER, DIMENSION(3) :: cell
389 38 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
390 38 : npgfb, nsgfa, nsgfb
391 38 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
392 : LOGICAL :: do_symmetric, found, trans
393 : REAL(KIND=dp) :: tab
394 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qab, tkab
395 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kab
396 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
397 38 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
398 38 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
399 342 : TYPE(block_p_type), DIMENSION(3, 2) :: blocks_tr
400 : TYPE(cell_type), POINTER :: qs_cell
401 38 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
402 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
403 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
404 :
405 : !$ INTEGER(kind=omp_lock_kind), &
406 38 : !$ ALLOCATABLE, DIMENSION(:) :: locks
407 : !$ INTEGER :: lock_num, hash, hash1, hash2
408 : !$ INTEGER(KIND=int_8) :: iatom8
409 : !$ INTEGER, PARAMETER :: nlock = 501
410 :
411 : MARK_USED(int_8)
412 :
413 38 : CALL timeset(routineN, handle)
414 :
415 38 : nkind = SIZE(qs_kind_set)
416 :
417 : ! check for symmetry
418 38 : CPASSERT(SIZE(sab_nl) > 0)
419 38 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
420 :
421 : ! prepare basis set
422 190 : ALLOCATE (basis_set_list(nkind))
423 38 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
424 :
425 : ! *** Allocate work storage ***
426 38 : ldsab = get_memory_usage(qs_kind_set, basis_type)
427 :
428 : CALL get_qs_env(qs_env=qs_env, &
429 : particle_set=particle_set, &
430 : cell=qs_cell, &
431 38 : natom=natom)
432 :
433 : !$OMP PARALLEL DEFAULT(NONE) &
434 : !$OMP SHARED (ldsab,do_symmetric, sab_nl, rc,&
435 : !$OMP ncoset,matrix_tr,matrix_rt,basis_set_list,qs_cell,natom,locks) &
436 : !$OMP PRIVATE (blocks_tr,kab,qab,tab,ikind,jkind,iatom,jatom,rab,rac,rbc,cell, &
437 : !$OMP basis_set_a, basis_set_b, nseta, ncoa, ncob, ltab, nsetb, tkab, &
438 : !$OMP irow, icol, found, trans, sgfa, sgfb, iset, jset, i, &
439 : !$OMP hash, hash1, hash2, iatom8, slot, lock_num) &
440 : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, zeta, scon_a) &
441 : !$OMP PRIVATE (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, zetb, scon_b) &
442 : !$OMP SHARED(particle_set) &
443 38 : !$OMP PRIVATE(ra, rb)
444 :
445 : !$OMP SINGLE
446 : !$ ALLOCATE (locks(nlock))
447 : !$OMP END SINGLE
448 :
449 : !$OMP DO
450 : !$ DO lock_num = 1, nlock
451 : !$ call omp_init_lock(locks(lock_num))
452 : !$ END DO
453 : !$OMP END DO
454 :
455 : ALLOCATE (kab(ldsab, ldsab, 3), qab(ldsab, ldsab))
456 :
457 : !$OMP DO SCHEDULE(GUIDED)
458 : DO slot = 1, sab_nl(1)%nl_size
459 :
460 : ikind = sab_nl(1)%nlist_task(slot)%ikind
461 : jkind = sab_nl(1)%nlist_task(slot)%jkind
462 : iatom = sab_nl(1)%nlist_task(slot)%iatom
463 : jatom = sab_nl(1)%nlist_task(slot)%jatom
464 : cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
465 : rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
466 :
467 : basis_set_a => basis_set_list(ikind)%gto_basis_set
468 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
469 : basis_set_b => basis_set_list(jkind)%gto_basis_set
470 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
471 :
472 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
473 : !$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
474 :
475 : ! basis ikind
476 : first_sgfa => basis_set_a%first_sgf
477 : la_max => basis_set_a%lmax
478 : la_min => basis_set_a%lmin
479 : npgfa => basis_set_a%npgf
480 : nsgfa => basis_set_a%nsgf_set
481 : rpgfa => basis_set_a%pgf_radius
482 : set_radius_a => basis_set_a%set_radius
483 : scon_a => basis_set_a%scon
484 : zeta => basis_set_a%zet
485 : ! basis jkind
486 : first_sgfb => basis_set_b%first_sgf
487 : lb_max => basis_set_b%lmax
488 : lb_min => basis_set_b%lmin
489 : npgfb => basis_set_b%npgf
490 : nsgfb => basis_set_b%nsgf_set
491 : rpgfb => basis_set_b%pgf_radius
492 : set_radius_b => basis_set_b%set_radius
493 : scon_b => basis_set_b%scon
494 : zetb => basis_set_b%zet
495 :
496 : nseta = basis_set_a%nset
497 : nsetb = basis_set_b%nset
498 :
499 : IF (do_symmetric) THEN
500 : IF (iatom <= jatom) THEN
501 : irow = iatom
502 : icol = jatom
503 : ELSE
504 : irow = jatom
505 : icol = iatom
506 : END IF
507 : ELSE
508 : irow = iatom
509 : icol = jatom
510 : END IF
511 : DO ir = 1, 3
512 : NULLIFY (blocks_tr(ir, 1)%block)
513 : CALL dbcsr_get_block_p(matrix=matrix_tr(ir)%matrix, &
514 : row=irow, col=icol, BLOCK=blocks_tr(ir, 1)%block, found=found)
515 : CPASSERT(found)
516 : NULLIFY (blocks_tr(ir, 2)%block)
517 : CALL dbcsr_get_block_p(matrix=matrix_rt(ir)%matrix, &
518 : row=irow, col=icol, BLOCK=blocks_tr(ir, 2)%block, found=found)
519 : CPASSERT(found)
520 : END DO
521 :
522 : ! The kinetic integrals depend only on rab (also for the screening)
523 : tab = NORM2(rab)
524 :
525 : ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
526 : ra = pbc(particle_set(iatom)%r(:), qs_cell)
527 : rb(:) = ra(:) + rab(:)
528 : rac = pbc(rc, ra, qs_cell)
529 : rbc = rac + rab
530 :
531 : trans = do_symmetric .AND. (iatom > jatom)
532 :
533 : DO iset = 1, nseta
534 :
535 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
536 : sgfa = first_sgfa(1, iset)
537 :
538 : DO jset = 1, nsetb
539 :
540 : IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE
541 :
542 : !$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
543 : !$ hash = MOD(hash1 + hash2, nlock) + 1
544 :
545 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
546 : sgfb = first_sgfb(1, jset)
547 :
548 : ! calculate integrals
549 : ltab = MAX(npgfa(iset)*ncoset(la_max(iset) + 1), npgfb(jset)*ncoset(lb_max(jset) + 1))
550 : ALLOCATE (tkab(ltab, ltab))
551 : CALL kinetic(la_max(iset) + 1, la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
552 : lb_max(jset) + 1, lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
553 : rab, tkab)
554 : ! Ordered products T*r and r*T
555 : DO i = 1, 2
556 : CALL apply_cartesian_operator(la_max(iset), la_min(iset), npgfa(iset), &
557 : lb_max(jset), lb_min(jset), npgfb(jset), tkab, kab, &
558 : right=(i == 1), derivative=.FALSE., origin=MERGE(rbc, rac, i == 1))
559 : DO ir = 1, 3
560 : CALL contraction(kab(:, :, ir), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
561 : cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), trans=trans)
562 :
563 : !$ CALL omp_set_lock(locks(hash))
564 : CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), blocks_tr(ir, i)%block, &
565 : sgfa, sgfb, trans=trans)
566 : !$ CALL omp_unset_lock(locks(hash))
567 : END DO
568 : END DO
569 : DEALLOCATE (tkab)
570 :
571 : END DO
572 : END DO
573 : END DO !iterator
574 : DEALLOCATE (kab, qab)
575 : !$OMP DO
576 : !$ DO lock_num = 1, nlock
577 : !$ call omp_destroy_lock(locks(lock_num))
578 : !$ END DO
579 : !$OMP END DO
580 :
581 : !$OMP SINGLE
582 : !$ DEALLOCATE (locks)
583 : !$OMP END SINGLE NOWAIT
584 :
585 : !$OMP END PARALLEL
586 :
587 : ! Release work storage
588 38 : DEALLOCATE (basis_set_list)
589 38 : CALL timestop(handle)
590 :
591 152 : END SUBROUTINE build_tr_matrices
592 :
593 : ! **************************************************************************************************
594 : !> \brief Build the ordered products V_core*r and r*V_core.
595 : !> \param matrix_rcore_hr ...
596 : !> \param matrix_rcore_rh ...
597 : !> \param qs_env ...
598 : !> \param qs_kind_set ...
599 : !> \param basis_type ...
600 : !> \param sab_nl ...
601 : !> \param rf ...
602 : !> \author Edward Ditler, Tomas Zimmermann
603 : ! **************************************************************************************************
604 38 : SUBROUTINE build_rcore_matrices(matrix_rcore_hr, matrix_rcore_rh, qs_env, qs_kind_set, basis_type, sab_nl, rf)
605 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_rcore_hr, matrix_rcore_rh
606 : TYPE(qs_environment_type), POINTER :: qs_env
607 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
608 : CHARACTER(LEN=*) :: basis_type
609 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
610 : POINTER :: sab_nl
611 : REAL(KIND=dp), DIMENSION(3) :: rf
612 :
613 : CHARACTER(len=*), PARAMETER :: routineN = 'build_rcore_matrices'
614 : INTEGER, PARAMETER :: nexp_max = 30
615 :
616 : INTEGER :: atom_a, atom_b, handle, i, iatom, icol, idir, ikind, inode, irow, iset, jatom, &
617 : jkind, jset, katom, kkind, ldai, ldsab, maxco, maxder, maxl, maxlgto, maxlppl, maxnset, &
618 : maxsgf, mepos, n_local, ncoa, ncob, nder, nexp_lpot, nexp_ppl, nimages, nkind, nloc, &
619 : nseta, nsetb, nthread, sgfa, sgfb
620 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
621 : INTEGER, DIMENSION(1:10) :: nrloc
622 : INTEGER, DIMENSION(3) :: cellind
623 38 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, &
624 38 : nct_lpot, npgfa, npgfb, nsgfa, nsgfb
625 38 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
626 : INTEGER, DIMENSION(nexp_max) :: nct_ppl
627 : LOGICAL :: do_symmetric, do_upf, dokp, ecp_local, &
628 : ecp_semi_local, found, lpotextended
629 : REAL(KIND=dp) :: alpha, dab, dac, dbc, ppl_radius
630 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: hab, qab
631 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ppl_work, rhab, work
632 : REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
633 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, raf, rb, rbc
634 : REAL(KIND=dp), DIMENSION(4, nexp_max) :: cval_ppl
635 76 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
636 38 : set_radius_a, set_radius_b
637 38 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, rpgfa, rpgfb, scon_a, scon_b, &
638 38 : sphi_a, sphi_b, zeta, zetb
639 : REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
640 38 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
641 342 : TYPE(block_p_type), DIMENSION(3, 2) :: blocks_rcore
642 : TYPE(cell_type), POINTER :: cell
643 : TYPE(gth_potential_type), POINTER :: gth_potential
644 38 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
645 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
646 : TYPE(neighbor_list_iterator_p_type), &
647 76 : DIMENSION(:), POINTER :: ap_iterator, nl_iterator
648 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
649 38 : POINTER :: sac_ppl
650 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
651 : TYPE(sgp_potential_type), POINTER :: sgp_potential
652 :
653 38 : CALL timeset(routineN, handle)
654 :
655 : CALL get_qs_env(qs_env=qs_env, &
656 : atomic_kind_set=atomic_kind_set, &
657 : qs_kind_set=qs_kind_set, &
658 : particle_set=particle_set, &
659 : sac_ppl=sac_ppl, &
660 38 : cell=cell)
661 :
662 : ! check for symmetry
663 38 : CPASSERT(SIZE(sab_nl) > 0)
664 38 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
665 :
666 38 : nkind = SIZE(qs_kind_set)
667 :
668 : ! prepare basis set
669 190 : ALLOCATE (basis_set_list(nkind))
670 38 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
671 :
672 38 : nder = 0
673 38 : nimages = 1
674 :
675 : alpha_ppl = 0
676 : nct_ppl = 0
677 : cval_ppl = 0
678 :
679 38 : nkind = SIZE(atomic_kind_set)
680 :
681 38 : dokp = (nimages > 1)
682 :
683 38 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
684 :
685 38 : maxder = ncoset(nder)
686 :
687 : CALL get_qs_kind_set(qs_kind_set, maxco=maxco, maxlgto=maxlgto, &
688 : maxsgf=maxsgf, maxnset=maxnset, maxlppl=maxlppl, &
689 38 : basis_type=basis_type)
690 :
691 38 : maxl = MAX(maxlgto, maxlppl)
692 38 : CALL init_orbital_pointers(2*maxl + 2*nder + 2)
693 :
694 : !tz: maxco in maxco*ncoset(maxlgto+1) is an overkill,
695 : ! properly there should be maxpgf*ncoset(maxlgto+1), but maxpgf is difficult to get
696 38 : ldsab = MAX(maxco, ncoset(maxlppl), maxsgf, maxlppl, maxco*ncoset(maxlgto + 1))
697 38 : ldai = ncoset(2*maxlgto + 2)
698 :
699 114 : DO ikind = 1, nkind
700 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type, &
701 76 : sgp_potential=sgp_potential)
702 76 : IF (ASSOCIATED(sgp_potential)) THEN
703 0 : IF (ALLOCATED(sgp_potential%upf_data)) THEN
704 0 : IF (sgp_potential%upf_local%radial%degree < 2*maxlgto + 2) THEN
705 : CALL upf_local_init(sgp_potential%upf_data, sgp_potential%alpha_core_charge, &
706 0 : 2*maxlgto + 2, sgp_potential%upf_local)
707 : END IF
708 : END IF
709 : END IF
710 114 : IF (ASSOCIATED(basis_set_a)) THEN
711 76 : basis_set_list(ikind)%gto_basis_set => basis_set_a
712 : ELSE
713 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
714 : END IF
715 : END DO
716 :
717 : nthread = 1
718 38 : !$ nthread = omp_get_max_threads()
719 :
720 38 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
721 :
722 : ! iterator for basis/potential list
723 38 : CALL neighbor_list_iterator_create(ap_iterator, sac_ppl, search=.TRUE., nthread=nthread)
724 :
725 : !$OMP PARALLEL &
726 : !$OMP DEFAULT (NONE) &
727 : !$OMP SHARED (nl_iterator, ap_iterator, basis_set_list, &
728 : !$OMP atomic_kind_set, qs_kind_set, particle_set, &
729 : !$OMP sab_nl, sac_ppl, nthread, ncoset, nkind, &
730 : !$OMP atom_of_kind, ldsab, maxnset, maxder, &
731 : !$OMP maxlgto, nder, maxco, dokp, cell) &
732 : !$OMP SHARED (matrix_rcore_hr, matrix_rcore_rh, rf) &
733 : !$OMP PRIVATE (ikind, jkind, inode, iatom, jatom, rab, basis_set_a, basis_set_b, atom_a) &
734 : !$OMP PRIVATE (atom_b) &
735 : !$OMP PRIVATE (nsetb) &
736 : !$OMP PRIVATE (dab, irow, icol, blocks_rcore, found, iset, ncoa, i) &
737 : !$OMP PRIVATE (sgfa, jset, ncob, sgfb, work, hab, rhab, kkind, nseta) &
738 : !$OMP PRIVATE (gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended) &
739 : !$OMP PRIVATE (ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl) &
740 : !$OMP PRIVATE (ecp_semi_local, nct_lpot, cval_ppl, cval_lpot, rac, dac, rbc, dbc) &
741 : !$OMP PRIVATE (mepos) &
742 : !$OMP PRIVATE (katom, ppl_work, cellind, ecp_local) &
743 : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, sphi_a, zeta, scon_a) &
744 : !$OMP PRIVATE (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, sphi_b, zetb, scon_b) &
745 : !$OMP PRIVATE (nloc, nrloc, aloc, bloc, n_local, a_local, c_local, ldai) &
746 38 : !$OMP PRIVATE (ra, rb, qab, raf, do_upf)
747 :
748 : mepos = 0
749 : !$ mepos = omp_get_thread_num()
750 :
751 : ALLOCATE (hab(ldsab, ldsab), rhab(ldsab, ldsab, 3), work(ldsab, ldsab*(nder + 1), 3))
752 : ALLOCATE (qab(ldsab, ldsab))
753 :
754 : ldai = ncoset(2*maxlgto + 2)
755 : ALLOCATE (ppl_work(ldai, ldai, MAX(maxder, 2*maxlgto + 2 + 1)))
756 :
757 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
758 :
759 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, inode=inode, &
760 : iatom=iatom, jatom=jatom, r=rab, cell=cellind)
761 :
762 : basis_set_a => basis_set_list(ikind)%gto_basis_set
763 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
764 : basis_set_b => basis_set_list(jkind)%gto_basis_set
765 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
766 :
767 : atom_a = atom_of_kind(iatom)
768 : atom_b = atom_of_kind(jatom)
769 :
770 : ! basis ikind
771 : first_sgfa => basis_set_a%first_sgf
772 : la_max => basis_set_a%lmax
773 : la_min => basis_set_a%lmin
774 : npgfa => basis_set_a%npgf
775 : nsgfa => basis_set_a%nsgf_set
776 : rpgfa => basis_set_a%pgf_radius
777 : set_radius_a => basis_set_a%set_radius
778 : sphi_a => basis_set_a%sphi
779 : zeta => basis_set_a%zet
780 : scon_a => basis_set_a%scon
781 : ! basis jkind
782 : first_sgfb => basis_set_b%first_sgf
783 : lb_max => basis_set_b%lmax
784 : lb_min => basis_set_b%lmin
785 : npgfb => basis_set_b%npgf
786 : nsgfb => basis_set_b%nsgf_set
787 : rpgfb => basis_set_b%pgf_radius
788 : set_radius_b => basis_set_b%set_radius
789 : sphi_b => basis_set_b%sphi
790 : zetb => basis_set_b%zet
791 : scon_b => basis_set_b%scon
792 :
793 : nseta = basis_set_a%nset
794 : nsetb = basis_set_b%nset
795 :
796 : ! *** Create matrix blocks for a new matrix block column ***
797 : irow = iatom
798 : icol = jatom
799 :
800 : DO idir = 1, 3
801 : NULLIFY (blocks_rcore(idir, 1)%block)
802 : CALL dbcsr_get_block_p(matrix=matrix_rcore_hr(idir)%matrix, &
803 : row=irow, col=icol, BLOCK=blocks_rcore(idir, 1)%block, found=found)
804 : CPASSERT(found)
805 : NULLIFY (blocks_rcore(idir, 2)%block)
806 : CALL dbcsr_get_block_p(matrix=matrix_rcore_rh(idir)%matrix, &
807 : row=irow, col=icol, BLOCK=blocks_rcore(idir, 2)%block, found=found)
808 : CPASSERT(found)
809 : END DO
810 :
811 : dab = NORM2(rab)
812 : ra = pbc(particle_set(iatom)%r(:), cell)
813 : rb(:) = ra(:) + rab(:)
814 :
815 : raf = pbc(rf, ra, cell)
816 :
817 : ! loop over all kinds for pseudopotential atoms
818 : DO kkind = 1, nkind
819 : CALL get_qs_kind(qs_kind_set(kkind), gth_potential=gth_potential, &
820 : sgp_potential=sgp_potential)
821 : do_upf = .FALSE.
822 : IF (ASSOCIATED(gth_potential)) THEN
823 : CALL get_potential(potential=gth_potential, &
824 : alpha_ppl=alpha, cexp_ppl=cexp_ppl, &
825 : lpot_present=lpotextended, ppl_radius=ppl_radius)
826 : nexp_ppl = 1
827 : alpha_ppl(1) = alpha
828 : nct_ppl(1) = SIZE(cexp_ppl)
829 : cval_ppl(1:nct_ppl(1), 1) = cexp_ppl(1:nct_ppl(1))
830 : IF (lpotextended) THEN
831 : CALL get_potential(potential=gth_potential, &
832 : nexp_lpot=nexp_lpot, alpha_lpot=alpha_lpot, nct_lpot=nct_lpot, cval_lpot=cval_lpot)
833 : CPASSERT(nexp_lpot < nexp_max)
834 : nexp_ppl = nexp_lpot + 1
835 : alpha_ppl(2:nexp_lpot + 1) = alpha_lpot(1:nexp_lpot)
836 : nct_ppl(2:nexp_lpot + 1) = nct_lpot(1:nexp_lpot)
837 : DO i = 1, nexp_lpot
838 : cval_ppl(1:nct_lpot(i), i + 1) = cval_lpot(1:nct_lpot(i), i)
839 : END DO
840 : END IF
841 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
842 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local, &
843 : ppl_radius=ppl_radius)
844 : do_upf = ALLOCATED(sgp_potential%upf_data)
845 : IF (do_upf) THEN
846 : CPASSERT(.NOT. ecp_local .AND. .NOT. ecp_semi_local)
847 : ELSE IF (ecp_local) THEN
848 : CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc)
849 : IF (SUM(ABS(aloc(1:nloc))) < 1.0e-12_dp) CYCLE
850 : nexp_ppl = nloc
851 : CPASSERT(nexp_ppl <= nexp_max)
852 : nct_ppl(1:nloc) = nrloc(1:nloc) - 1
853 : alpha_ppl(1:nloc) = bloc(1:nloc)
854 : cval_ppl(1, 1:nloc) = aloc(1:nloc)
855 : ELSE
856 : CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
857 : nexp_ppl = n_local
858 : CPASSERT(nexp_ppl <= nexp_max)
859 : nct_ppl(1:n_local) = 1
860 : alpha_ppl(1:n_local) = a_local(1:n_local)
861 : cval_ppl(1, 1:n_local) = c_local(1:n_local)
862 : END IF
863 : IF (ecp_semi_local) THEN
864 : CPABORT("VCD with semi-local ECPs not implemented")
865 : END IF
866 : ELSE
867 : CYCLE
868 : END IF
869 :
870 : CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
871 :
872 : DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
873 :
874 : CALL get_iterator_info(ap_iterator, mepos=mepos, jatom=katom, r=rac)
875 : dac = SQRT(SUM(rac*rac))
876 : rbc(:) = rac(:) - rab(:)
877 : dbc = SQRT(SUM(rbc*rbc))
878 : IF ((MAXVAL(set_radius_a(:)) + ppl_radius < dac) .OR. &
879 : (MAXVAL(set_radius_b(:)) + ppl_radius < dbc)) THEN
880 : CYCLE
881 : END IF
882 : DO iset = 1, nseta
883 : IF (set_radius_a(iset) + ppl_radius < dac) CYCLE
884 : ncoa = npgfa(iset)*ncoset(la_max(iset))
885 : ! ncoa = npgfa(iset)*(ncoset(la_max(iset))-ncoset(la_min(iset)-1))
886 : sgfa = first_sgfa(1, iset)
887 : DO jset = 1, nsetb
888 : IF (set_radius_b(jset) + ppl_radius < dbc) CYCLE
889 : ncob = npgfb(jset)*ncoset(lb_max(jset))
890 : ! ncob = npgfb(jset)*(ncoset(lb_max(jset))-ncoset(lb_min(jset)-1))
891 : sgfb = first_sgfb(1, jset)
892 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
893 : ! *** Calculate the GTH pseudo potential forces ***
894 : hab = 0
895 : rhab = 0
896 : ppl_work = 0
897 : work = 0
898 :
899 : IF (do_upf) THEN
900 : CALL upf_local_matrix(sgp_potential%upf_local, &
901 : la_max(iset) + 1, zeta(:npgfa(iset), iset), &
902 : lb_max(jset) + 1, zetb(:npgfb(jset), jset), rab, rac, &
903 : hab(:npgfa(iset)*ncoset(la_max(iset) + 1), :npgfb(jset)*ncoset(lb_max(jset) + 1)))
904 : ELSE
905 : CALL ppl_integral( &
906 : la_max(iset) + 1, la_min(iset), npgfa(iset), &
907 : rpgfa(:, iset), zeta(:, iset), &
908 : lb_max(jset) + 1, lb_min(jset), npgfb(jset), &
909 : rpgfb(:, jset), zetb(:, jset), &
910 : nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
911 : rab, dab, rac, dac, rbc, dbc, hab(:, :), ppl_work)
912 : END IF
913 :
914 : ! product with r
915 : CALL apply_cartesian_operator(la_max(iset), 0, npgfa(iset), lb_max(jset), 0, npgfb(jset), &
916 : hab, rhab, right=.FALSE., derivative=.FALSE., origin=raf)
917 :
918 : DO idir = 1, 3
919 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
920 : 1.0_dp, rhab(1, 1, idir), SIZE(rhab, 1), &
921 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
922 : 0.0_dp, work(1, 1, idir), SIZE(work, 1))
923 : !$OMP CRITICAL(h_block_critical)
924 : DO i = 1, 2
925 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
926 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
927 : work(1, 1, idir), SIZE(work, 1), &
928 : 1.0_dp, blocks_rcore(idir, i)%block(sgfa, sgfb), &
929 : SIZE(blocks_rcore(idir, i)%block, 1))
930 : END DO
931 : !$OMP END CRITICAL(h_block_critical)
932 : END DO
933 : END DO
934 : END DO
935 : END DO
936 : END DO
937 : END DO ! iterator
938 :
939 : DEALLOCATE (hab, rhab, work, ppl_work)
940 :
941 : !$OMP END PARALLEL
942 :
943 38 : CALL neighbor_list_iterator_release(ap_iterator)
944 38 : CALL neighbor_list_iterator_release(nl_iterator)
945 :
946 38 : DEALLOCATE (atom_of_kind, basis_set_list)
947 :
948 38 : CALL timestop(handle)
949 :
950 190 : END SUBROUTINE build_rcore_matrices
951 :
952 : ! **************************************************************************************************
953 : !> \brief Build the ordered Hartree+XC products V_HXC*r and r*V_HXC.
954 : !> \param matrix_hr ...
955 : !> \param matrix_rh ...
956 : !> \param qs_env ...
957 : !> \param rc ...
958 : !> \author Edward Ditler, Tomas Zimmermann
959 : ! **************************************************************************************************
960 38 : SUBROUTINE build_vhxc_matrices(matrix_hr, matrix_rh, qs_env, rc)
961 : TYPE(dbcsr_p_type), DIMENSION(:, :), &
962 : INTENT(INOUT), POINTER :: matrix_hr, matrix_rh
963 : TYPE(qs_environment_type), POINTER :: qs_env
964 : REAL(KIND=dp), DIMENSION(3) :: rc
965 :
966 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_vhxc_matrices'
967 : INTEGER, PARAMETER :: nspins = 1
968 :
969 : INTEGER :: handle, idir, ispin
970 : REAL(kind=dp) :: edisp
971 38 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
972 38 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_rvxc, matrix_rvxc_desymm
973 : TYPE(pw_env_type), POINTER :: pw_env
974 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
975 38 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace, v_tau_rspace
976 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
977 : TYPE(qs_energy_type), POINTER :: energy
978 : TYPE(qs_ks_env_type), POINTER :: ks_env
979 : TYPE(qs_rho_type), POINTER :: rho_struct
980 : TYPE(section_vals_type), POINTER :: input, xc_section
981 :
982 38 : CALL timeset(routineN, handle)
983 :
984 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks, &
985 : ks_env=ks_env, &
986 : pw_env=pw_env, &
987 : input=input, &
988 : v_hartree_rspace=v_hartree_rspace, &
989 38 : energy=energy)
990 :
991 38 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
992 :
993 38 : NULLIFY (matrix_rvxc, matrix_rvxc_desymm)
994 38 : CALL dbcsr_allocate_matrix_set(matrix_rvxc, nspins, 3)
995 38 : CALL dbcsr_allocate_matrix_set(matrix_rvxc_desymm, nspins, 3)
996 :
997 76 : DO ispin = 1, nspins
998 190 : DO idir = 1, 3
999 114 : CALL dbcsr_init_p(matrix_rvxc(ispin, idir)%matrix)
1000 114 : CALL dbcsr_init_p(matrix_rvxc_desymm(ispin, idir)%matrix)
1001 :
1002 114 : CALL dbcsr_copy(matrix_rvxc_desymm(ispin, idir)%matrix, matrix_hr(1, 1)%matrix)
1003 114 : CALL dbcsr_set(matrix_rvxc_desymm(ispin, idir)%matrix, 0._dp)
1004 :
1005 114 : CALL dbcsr_copy(matrix_rvxc(ispin, idir)%matrix, matrix_ks(ispin)%matrix)
1006 152 : CALL dbcsr_set(matrix_rvxc(ispin, idir)%matrix, 0.0_dp)
1007 : END DO
1008 : END DO
1009 :
1010 38 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1011 38 : CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
1012 :
1013 38 : NULLIFY (v_rspace)
1014 38 : NULLIFY (v_tau_rspace)
1015 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
1016 : vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=energy%exc, &
1017 : edisp=edisp, dispersion_env=qs_env%dispersion_env, &
1018 38 : just_energy=.FALSE.)
1019 :
1020 38 : IF (.NOT. ASSOCIATED(v_rspace)) THEN
1021 0 : ALLOCATE (v_rspace(nspins))
1022 0 : DO ispin = 1, nspins
1023 0 : CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
1024 0 : CALL pw_zero(v_rspace(ispin))
1025 : END DO
1026 : END IF
1027 :
1028 76 : DO ispin = 1, nspins
1029 38 : CALL pw_axpy(v_hartree_rspace, v_rspace(ispin), 1.0_dp/v_hartree_rspace%pw_grid%dvol)
1030 : CALL integrate_rv_rspace(v_rspace=v_rspace(ispin), hmat=matrix_rvxc(ispin, :), qs_env=qs_env, &
1031 38 : rc=rc)
1032 :
1033 190 : DO idir = 1, 3
1034 114 : CALL dbcsr_scale(matrix_rvxc(ispin, idir)%matrix, v_rspace(ispin)%pw_grid%dvol)
1035 114 : CALL dbcsr_desymmetrize(matrix_rvxc(ispin, idir)%matrix, matrix_rvxc_desymm(ispin, idir)%matrix)
1036 114 : CALL dbcsr_add(matrix_hr(ispin, idir)%matrix, matrix_rvxc_desymm(ispin, idir)%matrix, 1.0_dp, 1.0_dp)
1037 152 : CALL dbcsr_add(matrix_rh(ispin, idir)%matrix, matrix_rvxc_desymm(ispin, idir)%matrix, 1.0_dp, 1.0_dp)
1038 : END DO
1039 : END DO
1040 :
1041 : ! return pw grids
1042 76 : DO ispin = 1, nspins
1043 76 : CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1044 : END DO
1045 38 : DEALLOCATE (v_rspace)
1046 :
1047 38 : CALL dbcsr_deallocate_matrix_set(matrix_rvxc)
1048 38 : CALL dbcsr_deallocate_matrix_set(matrix_rvxc_desymm)
1049 :
1050 38 : CALL timestop(handle)
1051 :
1052 38 : END SUBROUTINE build_vhxc_matrices
1053 :
1054 : ! **************************************************************************************************
1055 : !> \brief Calculates the integrals < mu | r * V | nu >
1056 : !> There is no direction_Or argument, because the potentials commute with r
1057 : !> This routine uses integrate_pgf_product directly. It could probably be rewritten to use
1058 : !> the new task_list interface.
1059 : !> \param v_rspace ...
1060 : !> \param hmat ...
1061 : !> \param qs_env ...
1062 : !> \param rc ...
1063 : !> \author Edward Ditler
1064 : ! **************************************************************************************************
1065 38 : SUBROUTINE integrate_rv_rspace(v_rspace, hmat, qs_env, rc)
1066 :
1067 : TYPE(pw_r3d_rs_type) :: v_rspace
1068 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: hmat
1069 : TYPE(qs_environment_type), POINTER :: qs_env
1070 : REAL(KIND=dp), DIMENSION(3) :: rc
1071 :
1072 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_rv_rspace'
1073 :
1074 : CHARACTER(len=default_string_length) :: my_basis_type
1075 : INTEGER :: bcol, brow, handle, iatom, idir, igrid_level, ikind, ikind_old, ilevel, img, &
1076 : ipair, ipgf, ipgf_new, iset, iset_new, iset_old, itask, ithread, jatom, jkind, jkind_old, &
1077 : jpgf, jpgf_new, jset, jset_new, jset_old, ldsab, maxco, maxlgto, maxpgf, maxset, &
1078 : maxsgf_set, na1, na2, natom, nb1, nb2, ncoa, ncoa_full, ncob, ncob_full, nkind, nseta, &
1079 : nsetb, nthread, sgfa, sgfb
1080 38 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1081 38 : npgfb, nsgfa, nsgfb
1082 38 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1083 : LOGICAL :: atom_pair_changed, atom_pair_done, distributed_grids, found, has_threads, &
1084 : my_compute_tau, my_gapw, new_set_pair_coming
1085 : REAL(KIND=dp) :: eps_rho_rspace, f, prefactor, radius, &
1086 : zetp
1087 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rab_inv, rac, rb, rbc, rp
1088 38 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1089 38 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, sphi_a, sphi_b, work, zeta, zetb
1090 38 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: habt, rhab, workt
1091 38 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
1092 38 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1093 38 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: h_block
1094 : TYPE(cell_type), POINTER :: cell
1095 : TYPE(dbcsr_distribution_type) :: dist
1096 : TYPE(dft_control_type), POINTER :: dft_control
1097 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1098 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1099 : TYPE(mp_comm_type) :: group
1100 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1101 : TYPE(pw_env_type), POINTER :: pw_env
1102 38 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1103 38 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1104 : TYPE(task_list_type), POINTER :: task_list, task_list_soft
1105 38 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
1106 : TYPE(virial_type), POINTER :: virial
1107 :
1108 38 : CALL timeset(routineN, handle)
1109 :
1110 38 : my_compute_tau = .FALSE.
1111 38 : my_gapw = .FALSE.
1112 38 : my_basis_type = "ORB"
1113 :
1114 : ! get the task lists
1115 : CALL get_qs_env(qs_env=qs_env, &
1116 : task_list=task_list, &
1117 38 : task_list_soft=task_list_soft)
1118 38 : CPASSERT(ASSOCIATED(task_list))
1119 :
1120 : ! the information on the grids is provided through pw_env
1121 : ! pw_env has to be the parent env for the potential grid (input)
1122 : ! there is an option to provide an external grid
1123 38 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1124 38 : CPASSERT(ASSOCIATED(pw_env))
1125 :
1126 : ! get all the general information on the system we are working on
1127 : CALL get_qs_env(qs_env=qs_env, &
1128 : atomic_kind_set=atomic_kind_set, &
1129 : qs_kind_set=qs_kind_set, &
1130 : cell=cell, &
1131 : dft_control=dft_control, &
1132 : particle_set=particle_set, &
1133 : virial=virial, &
1134 38 : natom=natom)
1135 :
1136 38 : rab = 0._dp
1137 38 : rac = 0._dp
1138 : rbc = 0._dp
1139 :
1140 : ! short cuts to task list variables
1141 38 : tasks => task_list%tasks
1142 38 : atom_pair_send => task_list%atom_pair_send
1143 38 : atom_pair_recv => task_list%atom_pair_recv
1144 :
1145 38 : CPASSERT(ASSOCIATED(pw_env))
1146 38 : CALL pw_env_get(pw_env, rs_grids=rs_v)
1147 :
1148 : ! get mpi group from rs_v
1149 38 : group = rs_v(1)%desc%group
1150 :
1151 : ! assign from pw_env
1152 38 : gridlevel_info => pw_env%gridlevel_info
1153 :
1154 : ! transform the potential on the rs_multigrids
1155 38 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1156 :
1157 38 : nkind = SIZE(qs_kind_set)
1158 :
1159 : ! needs to be consistent with rho_rspace
1160 38 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1161 :
1162 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1163 : maxco=maxco, &
1164 : maxlgto=maxlgto, &
1165 : maxsgf_set=maxsgf_set, &
1166 38 : basis_type=my_basis_type)
1167 :
1168 38 : distributed_grids = .FALSE.
1169 190 : DO igrid_level = 1, gridlevel_info%ngrid_levels
1170 38 : IF (rs_v(igrid_level)%desc%distributed) THEN
1171 152 : distributed_grids = .TRUE.
1172 : END IF
1173 : END DO
1174 :
1175 : nthread = 1
1176 38 : !$ nthread = omp_get_max_threads()
1177 :
1178 : ! get maximum numbers
1179 38 : maxset = 0
1180 38 : maxpgf = 0
1181 114 : DO ikind = 1, nkind
1182 : CALL get_qs_kind(qs_kind_set(ikind), &
1183 76 : basis_set=orb_basis_set, basis_type=my_basis_type)
1184 :
1185 76 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
1186 :
1187 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1188 76 : npgf=npgfa, nset=nseta)
1189 :
1190 76 : maxset = MAX(nseta, maxset)
1191 190 : maxpgf = MAX(MAXVAL(npgfa), maxpgf)
1192 : END DO
1193 :
1194 38 : ldsab = MAX(maxco, maxsgf_set, maxpgf*ncoset(maxlgto + 1))
1195 :
1196 : ! *** Allocate work storage ***
1197 38 : NULLIFY (habt, workt)
1198 38 : CALL reallocate(habt, 1, ldsab, 1, ldsab, 0, nthread)
1199 38 : CALL reallocate(workt, 1, ldsab, 1, maxsgf_set, 0, nthread)
1200 190 : ALLOCATE (rhab(ldsab, ldsab, 3))
1201 :
1202 152 : ALLOCATE (h_block(3))
1203 :
1204 38 : ithread = 0
1205 38 : !$ ithread = omp_get_thread_num()
1206 38 : work => workt(:, :, ithread)
1207 38 : hab => habt(:, :, ithread)
1208 62358 : hab(:, :) = 0._dp
1209 :
1210 38 : iset_old = -1; jset_old = -1
1211 38 : ikind_old = -1; jkind_old = -1
1212 :
1213 : ! Here we loop over gridlevels first, finalising the matrix after each grid level is
1214 : ! completed. On each grid level, we loop over atom pairs, which will only access
1215 : ! a single block of each matrix, so with OpenMP, each matrix block is only touched
1216 : ! by a single thread for each grid level
1217 190 : loop_gridlevels: DO igrid_level = 1, gridlevel_info%ngrid_levels
1218 608 : DO idir = 1, 3
1219 456 : CALL dbcsr_work_create(hmat(idir)%matrix, work_mutable=.TRUE., n=nthread)
1220 456 : CALL dbcsr_get_info(hmat(idir)%matrix, distribution=dist)
1221 456 : CALL dbcsr_distribution_get(dist, has_threads=has_threads)
1222 456 : !$ IF (.NOT. has_threads) &
1223 608 : !$ CPABORT("No thread distribution defined.")
1224 : END DO
1225 :
1226 494 : loop_pairs: DO ipair = 1, task_list%npairs(igrid_level)
1227 2299 : loop_tasks: DO itask = task_list%taskstart(ipair, igrid_level), task_list%taskstop(ipair, igrid_level)
1228 1805 : ilevel = tasks(itask)%grid_level
1229 1805 : img = tasks(itask)%image
1230 1805 : iatom = tasks(itask)%iatom
1231 1805 : jatom = tasks(itask)%jatom
1232 1805 : iset = tasks(itask)%iset
1233 1805 : jset = tasks(itask)%jset
1234 1805 : ipgf = tasks(itask)%ipgf
1235 1805 : jpgf = tasks(itask)%jpgf
1236 1805 : CPASSERT(img == 1)
1237 :
1238 : ! At the start of a block of tasks, get atom data (and kind data, if needed)
1239 1805 : IF (itask == task_list%taskstart(ipair, igrid_level)) THEN
1240 :
1241 342 : ikind = particle_set(iatom)%atomic_kind%kind_number
1242 342 : jkind = particle_set(jatom)%atomic_kind%kind_number
1243 :
1244 342 : IF (iatom <= jatom) THEN
1245 228 : brow = iatom
1246 228 : bcol = jatom
1247 : ELSE
1248 114 : brow = jatom
1249 114 : bcol = iatom
1250 : END IF
1251 :
1252 342 : IF (ikind /= ikind_old) THEN
1253 : CALL get_qs_kind(qs_kind_set(ikind), &
1254 38 : basis_set=orb_basis_set, basis_type=my_basis_type)
1255 :
1256 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1257 : first_sgf=first_sgfa, &
1258 : lmax=la_max, &
1259 : lmin=la_min, &
1260 : npgf=npgfa, &
1261 : nset=nseta, &
1262 : nsgf_set=nsgfa, &
1263 : set_radius=set_radius_a, &
1264 : sphi=sphi_a, &
1265 38 : zet=zeta)
1266 : END IF
1267 :
1268 342 : IF (jkind /= jkind_old) THEN
1269 : CALL get_qs_kind(qs_kind_set(jkind), &
1270 228 : basis_set=orb_basis_set, basis_type=my_basis_type)
1271 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1272 : first_sgf=first_sgfb, &
1273 : lmax=lb_max, &
1274 : lmin=lb_min, &
1275 : npgf=npgfb, &
1276 : nset=nsetb, &
1277 : nsgf_set=nsgfb, &
1278 : set_radius=set_radius_b, &
1279 : sphi=sphi_b, &
1280 228 : zet=zetb)
1281 :
1282 : END IF
1283 :
1284 1368 : DO idir = 1, 3
1285 1026 : NULLIFY (h_block(idir)%block)
1286 1026 : CALL dbcsr_get_block_p(hmat(idir)%matrix, brow, bcol, h_block(idir)%block, found)
1287 1368 : CPASSERT(found)
1288 : END DO
1289 :
1290 : ikind_old = ikind
1291 : jkind_old = jkind
1292 :
1293 : atom_pair_changed = .TRUE.
1294 :
1295 : ELSE
1296 :
1297 : atom_pair_changed = .FALSE.
1298 :
1299 : END IF
1300 :
1301 1805 : IF (atom_pair_changed .OR. iset_old /= iset .OR. jset_old /= jset) THEN
1302 : ! We reuse the hab(:, :) array to put the new integrals in.
1303 :
1304 342 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1305 342 : ncoa_full = npgfa(iset)*ncoset(la_max(iset) + 1)
1306 342 : sgfa = first_sgfa(1, iset)
1307 342 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1308 342 : ncob_full = npgfb(jset)*ncoset(lb_max(jset) + 1)
1309 342 : sgfb = first_sgfb(1, jset)
1310 :
1311 342 : IF (iatom <= jatom) THEN
1312 162108 : hab(1:ncoa_full, 1:ncob_full) = 0._dp
1313 : ELSE
1314 53010 : hab(1:ncob_full, 1:ncoa_full) = 0._dp
1315 : END IF
1316 :
1317 : iset_old = iset
1318 : jset_old = jset
1319 :
1320 : END IF
1321 :
1322 7220 : rab = tasks(itask)%rab
1323 : ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
1324 1805 : ra = pbc(particle_set(iatom)%r(:), cell)
1325 7220 : rb(:) = ra(:) + rab(:)
1326 1805 : rac = pbc(rc, ra, cell)
1327 7220 : rbc = rac + rab
1328 :
1329 1805 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
1330 1805 : f = zetb(jpgf, jset)/zetp
1331 7220 : rp(:) = ra(:) + f*rab(:)
1332 :
1333 7220 : prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
1334 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1335 : lb_min=lb_min(jset), lb_max=lb_max(jset), &
1336 : ra=ra, rb=rb, rp=rp, &
1337 : zetp=zetp, eps=eps_rho_rspace, &
1338 1805 : prefactor=prefactor, cutoff=1.0_dp)
1339 :
1340 1805 : na1 = (ipgf - 1)*ncoset(la_max(iset) + 1) + 1
1341 1805 : na2 = ipgf*ncoset(la_max(iset) + 1)
1342 1805 : nb1 = (jpgf - 1)*ncoset(lb_max(jset) + 1) + 1
1343 1805 : nb2 = jpgf*ncoset(lb_max(jset) + 1)
1344 :
1345 1805 : IF (iatom <= jatom) THEN
1346 : CALL integrate_pgf_product( &
1347 : la_max(iset) + 1, zeta(ipgf, iset), la_min(iset), &
1348 : lb_max(jset) + 1, zetb(jpgf, jset), lb_min(jset), &
1349 : ra, rab, rs_v(igrid_level), &
1350 : hab, o1=na1 - 1, o2=nb1 - 1, &
1351 : radius=radius, &
1352 1216 : calculate_forces=.FALSE.)
1353 : ELSE
1354 2356 : rab_inv = -rab
1355 : CALL integrate_pgf_product( &
1356 : lb_max(jset) + 1, zetb(jpgf, jset), lb_min(jset), &
1357 : la_max(iset) + 1, zeta(ipgf, iset), la_min(iset), &
1358 : rb, rab_inv, rs_v(igrid_level), &
1359 : hab, o1=nb1 - 1, o2=na1 - 1, &
1360 : radius=radius, &
1361 589 : calculate_forces=.FALSE.)
1362 : END IF
1363 :
1364 1805 : new_set_pair_coming = .FALSE.
1365 1805 : atom_pair_done = .FALSE.
1366 1805 : IF (itask < task_list%taskstop(ipair, igrid_level)) THEN
1367 1463 : ilevel = tasks(itask + 1)%grid_level
1368 1463 : img = tasks(itask + 1)%image
1369 1463 : iatom = tasks(itask + 1)%iatom
1370 1463 : jatom = tasks(itask + 1)%jatom
1371 1463 : iset_new = tasks(itask + 1)%iset
1372 1463 : jset_new = tasks(itask + 1)%jset
1373 1463 : ipgf_new = tasks(itask + 1)%ipgf
1374 1463 : jpgf_new = tasks(itask + 1)%jpgf
1375 1463 : IF (iset_new /= iset .OR. jset_new /= jset) THEN
1376 : new_set_pair_coming = .TRUE.
1377 : END IF
1378 : ELSE
1379 : ! do not forget the last block
1380 : new_set_pair_coming = .TRUE.
1381 1805 : atom_pair_done = .TRUE.
1382 : END IF
1383 :
1384 342 : IF (new_set_pair_coming) THEN
1385 : ! Increase lx, ly, lz by one to account for the | r * b >
1386 342 : IF (iatom <= jatom) THEN
1387 : CALL apply_cartesian_operator(la_max(iset), 0, npgfa(iset), lb_max(jset), 0, npgfb(jset), &
1388 228 : hab, rhab, right=.FALSE., derivative=.FALSE., origin=rac)
1389 :
1390 : ELSE
1391 : CALL apply_cartesian_operator(lb_max(jset), 0, npgfb(jset), la_max(iset), 0, npgfa(iset), &
1392 114 : hab, rhab, right=.TRUE., derivative=.FALSE., origin=rac)
1393 : END IF
1394 :
1395 : ! contract the block into h if we're done with the current set pair
1396 1368 : DO idir = 1, 3
1397 1368 : IF (iatom <= jatom) THEN
1398 112860 : work = 0._dp
1399 301815 : work(1:ncoa, 1:nsgfb(jset)) = MATMUL(rhab(1:ncoa, 1:ncob, idir), sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
1400 : h_block(idir)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
1401 : h_block(idir)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
1402 69255 : MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
1403 : ELSE
1404 36252 : work(1:ncob, 1:nsgfa(iset)) = MATMUL(rhab(1:ncob, 1:ncoa, idir), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
1405 : h_block(idir)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
1406 : h_block(idir)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
1407 15903 : MATMUL(TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)), work(1:ncob, 1:nsgfa(iset)))
1408 : END IF
1409 : END DO
1410 : END IF
1411 :
1412 : END DO loop_tasks
1413 : END DO loop_pairs
1414 :
1415 646 : DO idir = 1, 3
1416 608 : CALL dbcsr_finalize(hmat(idir)%matrix)
1417 : END DO
1418 :
1419 : END DO loop_gridlevels
1420 :
1421 : ! *** Release work storage ***
1422 38 : DEALLOCATE (habt, rhab, workt, h_block)
1423 :
1424 38 : CALL timestop(handle)
1425 :
1426 76 : END SUBROUTINE integrate_rv_rspace
1427 :
1428 : ! **************************************************************************************************
1429 : !> \brief Builds the overlap derivative wrt nuclear velocities
1430 : !> dS/dV = < mu | r | nu > * (nu - mu)
1431 : !> \param qs_env ...
1432 : !> \param matrix_dsdv ...
1433 : !> \param deltaR ...
1434 : !> \param rcc ...
1435 : !> \author Edward Ditler
1436 : ! **************************************************************************************************
1437 6 : SUBROUTINE build_dSdV_matrix(qs_env, matrix_dsdv, deltaR, rcc)
1438 : TYPE(qs_environment_type), POINTER :: qs_env
1439 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_dsdv
1440 : REAL(KIND=dp), DIMENSION(:, :) :: deltaR
1441 : REAL(KIND=dp), DIMENSION(3) :: rcc
1442 :
1443 : CHARACTER(len=*), PARAMETER :: routineN = 'build_dSdV_matrix'
1444 :
1445 : INTEGER :: handle, i
1446 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, my_matrix_dsdv, &
1447 6 : my_matrix_dsdv2, my_matrix_mom
1448 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1449 6 : POINTER :: sab_all
1450 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1451 :
1452 6 : CALL timeset(routineN, handle)
1453 :
1454 6 : NULLIFY (my_matrix_mom, my_matrix_dsdv, my_matrix_dsdv2, sab_all, qs_kind_set)
1455 :
1456 : CALL get_qs_env(qs_env=qs_env, &
1457 : sab_all=sab_all, &
1458 : qs_kind_set=qs_kind_set, &
1459 6 : matrix_ks=matrix_ks)
1460 :
1461 : ! my_matrix_mom is needed because build_local_moment only works correctly with symmetric matrices
1462 6 : CALL dbcsr_allocate_matrix_set(my_matrix_dsdv, 3)
1463 6 : CALL dbcsr_allocate_matrix_set(my_matrix_dsdv2, 3)
1464 6 : CALL dbcsr_allocate_matrix_set(my_matrix_mom, 3)
1465 :
1466 24 : DO i = 1, 3
1467 18 : ALLOCATE (my_matrix_dsdv(i)%matrix)
1468 18 : ALLOCATE (my_matrix_dsdv2(i)%matrix)
1469 18 : ALLOCATE (my_matrix_mom(i)%matrix)
1470 :
1471 18 : CALL dbcsr_copy(my_matrix_dsdv(i)%matrix, matrix_dsdv(i)%matrix)
1472 18 : CALL dbcsr_copy(my_matrix_dsdv2(i)%matrix, matrix_dsdv(i)%matrix)
1473 18 : CALL dbcsr_copy(my_matrix_mom(i)%matrix, matrix_ks(1)%matrix)
1474 :
1475 18 : CALL dbcsr_set(my_matrix_dsdv2(i)%matrix, 0.0_dp)
1476 18 : CALL dbcsr_set(my_matrix_dsdv(i)%matrix, 0.0_dp)
1477 18 : CALL dbcsr_set(my_matrix_mom(i)%matrix, 0.0_dp)
1478 24 : CALL dbcsr_set(matrix_dsdv(i)%matrix, 0.0_dp)
1479 : END DO
1480 :
1481 : CALL build_local_moment_matrix(qs_env, my_matrix_mom, 1, ref_point=rcc, &
1482 6 : neighbor_image=.TRUE.)
1483 :
1484 24 : DO i = 1, 3
1485 18 : CALL dbcsr_desymmetrize(my_matrix_mom(i)%matrix, my_matrix_dsdv(i)%matrix)
1486 24 : CALL dbcsr_copy(my_matrix_dsdv2(i)%matrix, my_matrix_dsdv(i)%matrix)
1487 : END DO
1488 :
1489 : ! delta_nu^A <mu|r|nu>
1490 : CALL hr_mult_by_delta_3d(my_matrix_dsdv, qs_kind_set, "ORB", sab_all, &
1491 6 : deltaR, direction_Or=.TRUE.)
1492 : ! -delta_mu^A <mu|r|nu>
1493 : CALL hr_mult_by_delta_3d(my_matrix_dsdv2, qs_kind_set, "ORB", sab_all, &
1494 6 : deltaR, direction_Or=.FALSE.)
1495 24 : DO i = 1, 3
1496 18 : CALL dbcsr_copy(matrix_dsdv(i)%matrix, my_matrix_dsdv(i)%matrix)
1497 24 : CALL dbcsr_add(matrix_dsdv(i)%matrix, my_matrix_dsdv2(i)%matrix, 1.0_dp, -1.0_dp)
1498 : END DO
1499 :
1500 6 : CALL dbcsr_deallocate_matrix_set(my_matrix_dsdv)
1501 6 : CALL dbcsr_deallocate_matrix_set(my_matrix_dsdv2)
1502 6 : CALL dbcsr_deallocate_matrix_set(my_matrix_mom)
1503 :
1504 6 : CALL timestop(handle)
1505 :
1506 6 : END SUBROUTINE build_dSdV_matrix
1507 :
1508 : ! **************************************************************************************************
1509 : !> \brief Build the two ordered nonlocal second-position products.
1510 : !> \param matrix_rcomr ...
1511 : !> \param matrix_rrcom ...
1512 : !> \param qs_kind_set ...
1513 : !> \param sab_all ...
1514 : !> \param sap_ppnl ...
1515 : !> \param eps_ppnl ...
1516 : !> \param particle_set ...
1517 : !> \param cell ...
1518 : !> \author Edward Ditler, Tomas Zimmermann
1519 : ! **************************************************************************************************
1520 2 : SUBROUTINE build_com_rpnl_r(matrix_rcomr, matrix_rrcom, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, &
1521 : particle_set, cell)
1522 :
1523 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: matrix_rcomr, matrix_rrcom
1524 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1525 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1526 : POINTER :: sab_all, sap_ppnl
1527 : REAL(KIND=dp), INTENT(IN) :: eps_ppnl
1528 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1529 : POINTER :: particle_set
1530 : TYPE(cell_type), POINTER :: cell
1531 :
1532 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_com_rpnl_r'
1533 :
1534 : INTEGER :: handle, i, iab, iac, iatom, ibc, icol, &
1535 : ikind, irow, j, jatom, jkind, kac, &
1536 : kbc, kkind, na, natom, nb, nkind, np, &
1537 : slot
1538 : INTEGER, DIMENSION(3) :: cell_b
1539 : LOGICAL :: found, ppnl_present
1540 : REAL(KIND=dp), DIMENSION(3) :: rab
1541 2 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1542 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1543 54 : TYPE(block_p_type), DIMENSION(3, 3, 2) :: blocks_rvr
1544 2 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
1545 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1546 2 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1547 :
1548 : !$ INTEGER(kind=omp_lock_kind), &
1549 2 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1550 : !$ INTEGER :: lock_num, hash
1551 : !$ INTEGER, PARAMETER :: nlock = 501
1552 :
1553 2 : ppnl_present = ASSOCIATED(sap_ppnl)
1554 2 : IF (.NOT. ppnl_present) RETURN
1555 :
1556 2 : CALL timeset(routineN, handle)
1557 2 : nkind = SIZE(qs_kind_set)
1558 2 : natom = SIZE(particle_set)
1559 :
1560 : ! sap_int needs to be shared as multiple threads need to access this
1561 2 : NULLIFY (sap_int)
1562 14 : ALLOCATE (sap_int(nkind*nkind))
1563 10 : DO i = 1, nkind*nkind
1564 8 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1565 10 : sap_int(i)%nalist = 0
1566 : END DO
1567 :
1568 : ! We put zero as a reference point, because actually in the integrals we need two different ones:
1569 : ! < a | (r - R^\lambda_\beta) * [V, r_\alpha - R^\eta_\alpha] | b >
1570 : ! The first reference point can be added in a seperate step as the term will be
1571 : ! - R^\lambda_\beta * < a | [V, r_\alpha] | b >
1572 : ! = + R^\lambda_\beta * < a | [r_\alpha, V] | b >
1573 : ! The second reference point is not important, because it disappears in the commutator
1574 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder=2, moment_mode=.TRUE., &
1575 2 : particle_set=particle_set, cell=cell, refpoint=[0._dp, 0._dp, 0._dp])
1576 :
1577 : ! *** Set up a sorting index
1578 2 : CALL sap_sort(sap_int)
1579 :
1580 10 : ALLOCATE (basis_set(nkind))
1581 6 : DO ikind = 1, nkind
1582 4 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1583 6 : IF (ASSOCIATED(orb_basis_set)) THEN
1584 4 : basis_set(ikind)%gto_basis_set => orb_basis_set
1585 : ELSE
1586 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1587 : END IF
1588 : END DO
1589 :
1590 : ! *** All integrals needed have been calculated and stored in sap_int
1591 : ! *** We now calculate the ordered second-position product matrix elements
1592 :
1593 : !$OMP PARALLEL &
1594 : !$OMP DEFAULT (NONE) &
1595 : !$OMP SHARED (basis_set, matrix_rcomr, matrix_rrcom, &
1596 : !$OMP sap_int, nkind, eps_ppnl, locks, sab_all) &
1597 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1598 : !$OMP iab, irow, icol, blocks_rvr, lock_num, &
1599 : !$OMP found, iac, ibc, alist_ac, alist_bc, &
1600 : !$OMP na, np, nb, kkind, kac, kbc, i, j, &
1601 2 : !$OMP hash, natom, acint, bcint, achint, bchint)
1602 :
1603 : !$OMP SINGLE
1604 : !$ ALLOCATE (locks(nlock))
1605 : !$OMP END SINGLE
1606 :
1607 : !$OMP DO
1608 : !$ DO lock_num = 1, nlock
1609 : !$ call omp_init_lock(locks(lock_num))
1610 : !$ END DO
1611 : !$OMP END DO
1612 :
1613 : !$OMP DO SCHEDULE(GUIDED)
1614 :
1615 : DO slot = 1, sab_all(1)%nl_size
1616 :
1617 : ikind = sab_all(1)%nlist_task(slot)%ikind
1618 : jkind = sab_all(1)%nlist_task(slot)%jkind
1619 : iatom = sab_all(1)%nlist_task(slot)%iatom
1620 : jatom = sab_all(1)%nlist_task(slot)%jatom
1621 : cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1622 : rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1623 :
1624 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
1625 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
1626 : iab = ikind + nkind*(jkind - 1)
1627 :
1628 : ! *** Create matrix blocks for a new matrix block column ***
1629 : irow = iatom
1630 : icol = jatom
1631 : DO i = 1, 3
1632 : DO j = 1, 3
1633 : ! t_alpha = MOD(i - 1, 3) + 1
1634 : ! t_beta = FLOOR(REAL(i - 1, dp)/3._dp) + 1
1635 :
1636 : CALL dbcsr_get_block_p(matrix_rcomr(i, j)%matrix, irow, icol, blocks_rvr(i, j, 1)%block, found)
1637 : CPASSERT(found)
1638 : CALL dbcsr_get_block_p(matrix_rrcom(i, j)%matrix, irow, icol, blocks_rvr(i, j, 2)%block, found)
1639 : CPASSERT(found)
1640 : END DO
1641 : END DO
1642 :
1643 : ! loop over all kinds for projector atom
1644 : DO kkind = 1, nkind
1645 : iac = ikind + nkind*(kkind - 1)
1646 : ibc = jkind + nkind*(kkind - 1)
1647 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
1648 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
1649 : CALL get_alist(sap_int(iac), alist_ac, iatom)
1650 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
1651 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
1652 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
1653 : DO kac = 1, alist_ac%nclist
1654 : DO kbc = 1, alist_bc%nclist
1655 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
1656 :
1657 : ! Some documention, considering the C1 O1 O2 molecule
1658 : ! The integrals are <a|p>
1659 : ! sap_int(1:2, 1:2) -> [(a=C, p=C), (a=C, p=O), (a=O, p=C), (a=O, p=0)]
1660 : ! the (a=O, p=C) entry has an alist
1661 : ! alist has two elements: O1, O2
1662 : ! alist(O1) -> clist(C1)%acint are the integrals
1663 : ! alist(O2) -> clist(C1)%acint are the integrals
1664 :
1665 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1666 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
1667 : acint => alist_ac%clist(kac)%acint
1668 : bcint => alist_bc%clist(kbc)%acint
1669 : achint => alist_ac%clist(kac)%achint
1670 : bchint => alist_bc%clist(kbc)%achint
1671 : na = SIZE(acint, 1)
1672 : np = SIZE(acint, 2)
1673 : nb = SIZE(bcint, 1)
1674 : !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1675 : !$ CALL omp_set_lock(locks(hash))
1676 :
1677 : DO i = 1, 3
1678 : DO j = 1, 3
1679 : ! matrix_rcomr(i, j) = (Vnl*r_i)*r_j - (r_i*Vnl)*r_j
1680 : blocks_rvr(i, j, 1)%block(1:na, 1:nb) = blocks_rvr(i, j, 1)%block(1:na, 1:nb) + &
1681 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_rr(i, j)))) - &
1682 : MATMUL(achint(1:na, 1:np, bi_r(i)), TRANSPOSE(bcint(1:nb, 1:np, bi_r(j))))
1683 :
1684 : ! matrix_rrcom(i, j) = r_j*(Vnl*r_i) - r_j*(r_i*Vnl)
1685 : blocks_rvr(i, j, 2)%block(1:na, 1:nb) = blocks_rvr(i, j, 2)%block(1:na, 1:nb) + &
1686 : MATMUL(achint(1:na, 1:np, bi_r(j)), TRANSPOSE(bcint(1:nb, 1:np, bi_r(i)))) - &
1687 : MATMUL(achint(1:na, 1:np, bi_rr(i, j)), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1688 : END DO
1689 : END DO
1690 :
1691 : !$ CALL omp_unset_lock(locks(hash))
1692 : EXIT ! We have found a match and there can be only one single match
1693 : END IF
1694 : END DO
1695 : END DO
1696 : END DO
1697 : DO i = 1, 3
1698 : DO j = 1, 3
1699 : NULLIFY (blocks_rvr(i, j, 1)%block)
1700 : NULLIFY (blocks_rvr(i, j, 2)%block)
1701 : END DO
1702 : END DO
1703 : END DO
1704 :
1705 : !$OMP DO
1706 : !$ DO lock_num = 1, nlock
1707 : !$ call omp_destroy_lock(locks(lock_num))
1708 : !$ END DO
1709 : !$OMP END DO
1710 :
1711 : !$OMP SINGLE
1712 : !$ DEALLOCATE (locks)
1713 : !$OMP END SINGLE NOWAIT
1714 :
1715 : !$OMP END PARALLEL
1716 :
1717 2 : CALL release_sap_int(sap_int)
1718 :
1719 2 : DEALLOCATE (basis_set)
1720 :
1721 2 : CALL timestop(handle)
1722 :
1723 6 : END SUBROUTINE build_com_rpnl_r
1724 :
1725 : ! **************************************************************************************************
1726 : !> \brief Calculate the double commutator [[Vnl, r], r]
1727 : !> \param matrix_rv ...
1728 : !> \param qs_kind_set ...
1729 : !> \param sab_orb ...
1730 : !> \param sap_ppnl ...
1731 : !> \param eps_ppnl ...
1732 : !> \param particle_set ...
1733 : !> \param pseudoatom Only consider pseudopotentials on atom lambda
1734 : !> \author Edward Ditler
1735 : ! **************************************************************************************************
1736 6 : SUBROUTINE build_dcom_rpnl(matrix_rv, qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, pseudoatom)
1737 :
1738 : TYPE(dbcsr_p_type), DIMENSION(:, :) :: matrix_rv
1739 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1740 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1741 : POINTER :: sab_orb, sap_ppnl
1742 : REAL(KIND=dp), INTENT(IN) :: eps_ppnl
1743 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1744 : POINTER :: particle_set
1745 : INTEGER, INTENT(IN) :: pseudoatom
1746 :
1747 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dcom_rpnl'
1748 :
1749 : INTEGER :: handle, i, iab, iac, iatom, ibc, icol, &
1750 : ikind, irow, j, jatom, jkind, kac, &
1751 : kbc, kkind, na, natom, nb, nkind, np, &
1752 : slot
1753 : INTEGER, DIMENSION(3) :: cell_b
1754 : LOGICAL :: found, ppnl_present
1755 : REAL(KIND=dp), DIMENSION(3) :: rab
1756 6 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1757 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1758 78 : TYPE(block_p_type), DIMENSION(3, 3) :: blocks_rv
1759 6 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
1760 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1761 6 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1762 :
1763 : !$ INTEGER(kind=omp_lock_kind), &
1764 6 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1765 : !$ INTEGER :: lock_num, hash
1766 : !$ INTEGER, PARAMETER :: nlock = 501
1767 :
1768 6 : ppnl_present = ASSOCIATED(sap_ppnl)
1769 6 : IF (.NOT. ppnl_present) RETURN
1770 :
1771 6 : CALL timeset(routineN, handle)
1772 6 : nkind = SIZE(qs_kind_set)
1773 6 : natom = SIZE(particle_set)
1774 :
1775 : ! sap_int needs to be shared as multiple threads need to access this
1776 6 : NULLIFY (sap_int)
1777 42 : ALLOCATE (sap_int(nkind*nkind))
1778 30 : DO i = 1, nkind*nkind
1779 24 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1780 30 : sap_int(i)%nalist = 0
1781 : END DO
1782 :
1783 : ! "nder" in moment_mode is "order"
1784 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder=2, moment_mode=.TRUE., &
1785 6 : particle_set=particle_set, pseudoatom=pseudoatom)
1786 :
1787 : ! *** Set up a sorting index
1788 6 : CALL sap_sort(sap_int)
1789 :
1790 30 : ALLOCATE (basis_set(nkind))
1791 18 : DO ikind = 1, nkind
1792 12 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1793 18 : IF (ASSOCIATED(orb_basis_set)) THEN
1794 12 : basis_set(ikind)%gto_basis_set => orb_basis_set
1795 : ELSE
1796 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1797 : END IF
1798 : END DO
1799 :
1800 : ! *** All integrals needed have been calculated and stored in sap_int
1801 : ! *** We now calculate the commutator matrix elements
1802 :
1803 : !$OMP PARALLEL &
1804 : !$OMP DEFAULT (NONE) &
1805 : !$OMP SHARED (basis_set, matrix_rv, &
1806 : !$OMP sap_int, nkind, eps_ppnl, locks, sab_orb) &
1807 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1808 : !$OMP iab, irow, icol, blocks_rv, lock_num, &
1809 : !$OMP found, iac, ibc, alist_ac, alist_bc, &
1810 : !$OMP na, np, nb, kkind, kac, kbc, i, j, &
1811 6 : !$OMP hash, natom, acint, bcint, achint, bchint)
1812 :
1813 : !$OMP SINGLE
1814 : !$ ALLOCATE (locks(nlock))
1815 : !$OMP END SINGLE
1816 :
1817 : !$OMP DO
1818 : !$ DO lock_num = 1, nlock
1819 : !$ call omp_init_lock(locks(lock_num))
1820 : !$ END DO
1821 : !$OMP END DO
1822 :
1823 : !$OMP DO SCHEDULE(GUIDED)
1824 :
1825 : DO slot = 1, sab_orb(1)%nl_size
1826 :
1827 : ikind = sab_orb(1)%nlist_task(slot)%ikind
1828 : jkind = sab_orb(1)%nlist_task(slot)%jkind
1829 : iatom = sab_orb(1)%nlist_task(slot)%iatom
1830 : jatom = sab_orb(1)%nlist_task(slot)%jatom
1831 : cell_b(:) = sab_orb(1)%nlist_task(slot)%cell(:)
1832 : rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
1833 :
1834 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
1835 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
1836 : iab = ikind + nkind*(jkind - 1)
1837 :
1838 : ! *** Create matrix blocks for a new matrix block column ***
1839 : IF (iatom <= jatom) THEN
1840 : irow = iatom
1841 : icol = jatom
1842 : ELSE
1843 : irow = jatom
1844 : icol = iatom
1845 : END IF
1846 :
1847 : DO i = 1, 3
1848 : DO j = 1, 3
1849 : CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, blocks_rv(i, j)%block, found)
1850 : blocks_rv(i, j)%block = 0._dp
1851 : CPASSERT(found)
1852 : END DO
1853 : END DO
1854 :
1855 : ! loop over all kinds for projector atom
1856 : DO kkind = 1, nkind
1857 : iac = ikind + nkind*(kkind - 1)
1858 : ibc = jkind + nkind*(kkind - 1)
1859 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
1860 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
1861 : CALL get_alist(sap_int(iac), alist_ac, iatom)
1862 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
1863 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
1864 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
1865 : DO kac = 1, alist_ac%nclist
1866 : DO kbc = 1, alist_bc%nclist
1867 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
1868 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1869 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
1870 : acint => alist_ac%clist(kac)%acint
1871 : bcint => alist_bc%clist(kbc)%acint
1872 : achint => alist_ac%clist(kac)%achint
1873 : bchint => alist_bc%clist(kbc)%achint
1874 : na = SIZE(acint, 1)
1875 : np = SIZE(acint, 2)
1876 : nb = SIZE(bcint, 1)
1877 : !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1878 : !$ CALL omp_set_lock(locks(hash))
1879 : ! Template:
1880 : ! blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
1881 : ! MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xV
1882 : IF (iatom <= jatom) THEN
1883 : ! r_alpha*Vnl*r_beta
1884 : blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) + &
1885 : MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
1886 :
1887 : blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) + &
1888 : MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
1889 :
1890 : blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) + &
1891 : MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
1892 :
1893 : blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) + &
1894 : MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
1895 :
1896 : blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) + &
1897 : MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
1898 :
1899 : blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) + &
1900 : MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
1901 :
1902 : blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) + &
1903 : MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
1904 :
1905 : blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) + &
1906 : MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
1907 :
1908 : blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) + &
1909 : MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
1910 :
1911 : ! -r_alpha*r_beta*Vnl
1912 : blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) - &
1913 : MATMUL(achint(1:na, 1:np, bi_xx), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1914 :
1915 : blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) - &
1916 : MATMUL(achint(1:na, 1:np, bi_xy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1917 :
1918 : blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) - &
1919 : MATMUL(achint(1:na, 1:np, bi_xz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1920 :
1921 : blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) - &
1922 : MATMUL(achint(1:na, 1:np, bi_xy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1923 :
1924 : blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) - &
1925 : MATMUL(achint(1:na, 1:np, bi_yy), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1926 :
1927 : blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) - &
1928 : MATMUL(achint(1:na, 1:np, bi_yz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1929 :
1930 : blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) - &
1931 : MATMUL(achint(1:na, 1:np, bi_xz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1932 :
1933 : blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) - &
1934 : MATMUL(achint(1:na, 1:np, bi_yz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1935 :
1936 : blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) - &
1937 : MATMUL(achint(1:na, 1:np, bi_zz), TRANSPOSE(bcint(1:nb, 1:np, bi_1)))
1938 :
1939 : ! -Vnl*r_beta*r_alpha
1940 : blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) - &
1941 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xx)))
1942 :
1943 : blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) - &
1944 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xy)))
1945 :
1946 : blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) - &
1947 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xz)))
1948 :
1949 : blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) - &
1950 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xy)))
1951 :
1952 : blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) - &
1953 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yy)))
1954 :
1955 : blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) - &
1956 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yz)))
1957 :
1958 : blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) - &
1959 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_xz)))
1960 :
1961 : blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) - &
1962 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_yz)))
1963 :
1964 : blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) - &
1965 : MATMUL(achint(1:na, 1:np, bi_1), TRANSPOSE(bcint(1:nb, 1:np, bi_zz)))
1966 :
1967 : ! +r_beta*Vnl*r_alpha
1968 : blocks_rv(1, 1)%block(1:na, 1:nb) = blocks_rv(1, 1)%block(1:na, 1:nb) + &
1969 : MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
1970 :
1971 : blocks_rv(1, 2)%block(1:na, 1:nb) = blocks_rv(1, 2)%block(1:na, 1:nb) + &
1972 : MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
1973 :
1974 : blocks_rv(1, 3)%block(1:na, 1:nb) = blocks_rv(1, 3)%block(1:na, 1:nb) + &
1975 : MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_x)))
1976 :
1977 : blocks_rv(2, 1)%block(1:na, 1:nb) = blocks_rv(2, 1)%block(1:na, 1:nb) + &
1978 : MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
1979 :
1980 : blocks_rv(2, 2)%block(1:na, 1:nb) = blocks_rv(2, 2)%block(1:na, 1:nb) + &
1981 : MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
1982 :
1983 : blocks_rv(2, 3)%block(1:na, 1:nb) = blocks_rv(2, 3)%block(1:na, 1:nb) + &
1984 : MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_y)))
1985 :
1986 : blocks_rv(3, 1)%block(1:na, 1:nb) = blocks_rv(3, 1)%block(1:na, 1:nb) + &
1987 : MATMUL(achint(1:na, 1:np, bi_x), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
1988 :
1989 : blocks_rv(3, 2)%block(1:na, 1:nb) = blocks_rv(3, 2)%block(1:na, 1:nb) + &
1990 : MATMUL(achint(1:na, 1:np, bi_y), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
1991 :
1992 : blocks_rv(3, 3)%block(1:na, 1:nb) = blocks_rv(3, 3)%block(1:na, 1:nb) + &
1993 : MATMUL(achint(1:na, 1:np, bi_z), TRANSPOSE(bcint(1:nb, 1:np, bi_z)))
1994 : ELSE
1995 : ! r_alpha*Vnl*r_beta
1996 : blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) + &
1997 : MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_x)))
1998 :
1999 : blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) + &
2000 : MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_y)))
2001 :
2002 : blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) + &
2003 : MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_z)))
2004 :
2005 : blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) + &
2006 : MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_x)))
2007 :
2008 : blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) + &
2009 : MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_y)))
2010 :
2011 : blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) + &
2012 : MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_z)))
2013 :
2014 : blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) + &
2015 : MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_x)))
2016 :
2017 : blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) + &
2018 : MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_y)))
2019 :
2020 : blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) + &
2021 : MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_z)))
2022 :
2023 : ! -r_alpha*r_beta*Vnl
2024 : blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) - &
2025 : MATMUL(bchint(1:nb, 1:np, bi_xx), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2026 :
2027 : blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) - &
2028 : MATMUL(bchint(1:nb, 1:np, bi_xy), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2029 :
2030 : blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) - &
2031 : MATMUL(bchint(1:nb, 1:np, bi_xz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2032 :
2033 : blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) - &
2034 : MATMUL(bchint(1:nb, 1:np, bi_xy), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2035 :
2036 : blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) - &
2037 : MATMUL(bchint(1:nb, 1:np, bi_yy), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2038 :
2039 : blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) - &
2040 : MATMUL(bchint(1:nb, 1:np, bi_yz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2041 :
2042 : blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) - &
2043 : MATMUL(bchint(1:nb, 1:np, bi_xz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2044 :
2045 : blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) - &
2046 : MATMUL(bchint(1:nb, 1:np, bi_yz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2047 :
2048 : blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) - &
2049 : MATMUL(bchint(1:nb, 1:np, bi_zz), TRANSPOSE(acint(1:na, 1:np, bi_1)))
2050 :
2051 : ! -Vnl*r_beta*r_alpha
2052 : blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) - &
2053 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xx)))
2054 :
2055 : blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) - &
2056 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xy)))
2057 :
2058 : blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) - &
2059 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xz)))
2060 :
2061 : blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) - &
2062 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xy)))
2063 :
2064 : blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) - &
2065 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_yy)))
2066 :
2067 : blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) - &
2068 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_yz)))
2069 :
2070 : blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) - &
2071 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_xz)))
2072 :
2073 : blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) - &
2074 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_yz)))
2075 :
2076 : blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) - &
2077 : MATMUL(bchint(1:nb, 1:np, bi_1), TRANSPOSE(acint(1:na, 1:np, bi_zz)))
2078 :
2079 : ! +r_beta*Vnl*r_alpha
2080 : blocks_rv(1, 1)%block(1:nb, 1:na) = blocks_rv(1, 1)%block(1:nb, 1:na) + &
2081 : MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_x)))
2082 :
2083 : blocks_rv(1, 2)%block(1:nb, 1:na) = blocks_rv(1, 2)%block(1:nb, 1:na) + &
2084 : MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_x)))
2085 :
2086 : blocks_rv(1, 3)%block(1:nb, 1:na) = blocks_rv(1, 3)%block(1:nb, 1:na) + &
2087 : MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_x)))
2088 :
2089 : blocks_rv(2, 1)%block(1:nb, 1:na) = blocks_rv(2, 1)%block(1:nb, 1:na) + &
2090 : MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_y)))
2091 :
2092 : blocks_rv(2, 2)%block(1:nb, 1:na) = blocks_rv(2, 2)%block(1:nb, 1:na) + &
2093 : MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_y)))
2094 :
2095 : blocks_rv(2, 3)%block(1:nb, 1:na) = blocks_rv(2, 3)%block(1:nb, 1:na) + &
2096 : MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_y)))
2097 :
2098 : blocks_rv(3, 1)%block(1:nb, 1:na) = blocks_rv(3, 1)%block(1:nb, 1:na) + &
2099 : MATMUL(bchint(1:nb, 1:np, bi_x), TRANSPOSE(acint(1:na, 1:np, bi_z)))
2100 :
2101 : blocks_rv(3, 2)%block(1:nb, 1:na) = blocks_rv(3, 2)%block(1:nb, 1:na) + &
2102 : MATMUL(bchint(1:nb, 1:np, bi_y), TRANSPOSE(acint(1:na, 1:np, bi_z)))
2103 :
2104 : blocks_rv(3, 3)%block(1:nb, 1:na) = blocks_rv(3, 3)%block(1:nb, 1:na) + &
2105 : MATMUL(bchint(1:nb, 1:np, bi_z), TRANSPOSE(acint(1:na, 1:np, bi_z)))
2106 :
2107 : END IF
2108 : !$ CALL omp_unset_lock(locks(hash))
2109 : EXIT ! We have found a match and there can be only one single match
2110 : END IF
2111 : END DO
2112 : END DO
2113 : END DO
2114 : DO i = 1, 3
2115 : NULLIFY (blocks_rv(i, 1)%block)
2116 : NULLIFY (blocks_rv(i, 2)%block)
2117 : NULLIFY (blocks_rv(i, 3)%block)
2118 : END DO
2119 : END DO
2120 :
2121 : !$OMP DO
2122 : !$ DO lock_num = 1, nlock
2123 : !$ call omp_destroy_lock(locks(lock_num))
2124 : !$ END DO
2125 : !$OMP END DO
2126 :
2127 : !$OMP SINGLE
2128 : !$ DEALLOCATE (locks)
2129 : !$OMP END SINGLE NOWAIT
2130 :
2131 : !$OMP END PARALLEL
2132 :
2133 6 : CALL release_sap_int(sap_int)
2134 :
2135 6 : DEALLOCATE (basis_set)
2136 :
2137 6 : CALL timestop(handle)
2138 :
2139 18 : END SUBROUTINE build_dcom_rpnl
2140 :
2141 : ! **************************************************************************************************
2142 : !> \brief Apply the operator \delta_\mu^\lambda to zero out all elements of the matrix
2143 : !> which don't fulfill the condition.
2144 : !> Operates on matrix_hr(1:3) instead of a single matrix
2145 : !> \param matrix_hr ...
2146 : !> \param qs_kind_set ...
2147 : !> \param basis_type ...
2148 : !> \param sab_nl ...
2149 : !> \param deltaR ...
2150 : !> \param direction_Or ...
2151 : !> \author Edward Ditler
2152 : ! **************************************************************************************************
2153 18 : SUBROUTINE hr_mult_by_delta_3d(matrix_hr, qs_kind_set, basis_type, sab_nl, deltaR, direction_Or)
2154 :
2155 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hr
2156 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2157 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
2158 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2159 : POINTER :: sab_nl
2160 : REAL(KIND=dp), DIMENSION(:, :) :: deltaR
2161 : LOGICAL :: direction_Or
2162 :
2163 : CHARACTER(len=*), PARAMETER :: routineN = 'hr_mult_by_delta_3d'
2164 :
2165 : INTEGER :: handle, iatom, icol, ikind, ir, irow, &
2166 : jatom, jkind, ldsab, mepos, nkind, &
2167 : nseta, nsetb, nthread
2168 : INTEGER, DIMENSION(3) :: cell
2169 18 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
2170 18 : npgfb, nsgfa, nsgfb
2171 18 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
2172 : LOGICAL :: do_symmetric, found
2173 : REAL(KIND=dp), DIMENSION(3) :: rab
2174 18 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
2175 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: kx_block, ky_block, kz_block, rpgfa, &
2176 18 : rpgfb, scon_a, scon_b, sphi_a, sphi_b, &
2177 18 : zeta, zetb
2178 18 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
2179 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
2180 : TYPE(neighbor_list_iterator_p_type), &
2181 18 : DIMENSION(:), POINTER :: nl_iterator
2182 :
2183 18 : CALL timeset(routineN, handle)
2184 :
2185 18 : nkind = SIZE(qs_kind_set)
2186 :
2187 : ! check for symmetry
2188 18 : CPASSERT(SIZE(sab_nl) > 0)
2189 18 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
2190 :
2191 : ! prepare basis set
2192 90 : ALLOCATE (basis_set_list(nkind))
2193 18 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
2194 :
2195 : ! *** Allocate work storage ***
2196 18 : ldsab = get_memory_usage(qs_kind_set, basis_type)
2197 :
2198 : nthread = 1
2199 18 : !$ nthread = omp_get_max_threads()
2200 : ! Iterate of neighbor list
2201 18 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
2202 :
2203 : !$OMP PARALLEL DEFAULT(NONE) &
2204 : !$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric) &
2205 : !$OMP SHARED (ncoset,matrix_hr,basis_set_list) &
2206 : !$OMP SHARED (direction_or, deltar) &
2207 : !$OMP PRIVATE (kx_block,ky_block,kz_block,mepos,ikind,jkind,iatom,jatom,rab,cell) &
2208 : !$OMP PRIVATE (basis_set_a,basis_set_b) &
2209 : !$OMP PRIVATE (nseta, nsetb) &
2210 : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, rpgfa, set_radius_a, sphi_a, zeta, scon_a) &
2211 : !$OMP PRIVATE (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, set_radius_b, sphi_b, zetb, scon_b) &
2212 18 : !$OMP PRIVATE (irow, icol, found)
2213 :
2214 : mepos = 0
2215 : !$ mepos = omp_get_thread_num()
2216 :
2217 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
2218 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
2219 : iatom=iatom, jatom=jatom, r=rab, cell=cell)
2220 : basis_set_a => basis_set_list(ikind)%gto_basis_set
2221 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
2222 : basis_set_b => basis_set_list(jkind)%gto_basis_set
2223 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
2224 : ! basis ikind
2225 : first_sgfa => basis_set_a%first_sgf
2226 : la_max => basis_set_a%lmax
2227 : la_min => basis_set_a%lmin
2228 : npgfa => basis_set_a%npgf
2229 : nsgfa => basis_set_a%nsgf_set
2230 : rpgfa => basis_set_a%pgf_radius
2231 : set_radius_a => basis_set_a%set_radius
2232 : sphi_a => basis_set_a%sphi
2233 : zeta => basis_set_a%zet
2234 : scon_a => basis_set_a%scon
2235 : ! basis jkind
2236 : first_sgfb => basis_set_b%first_sgf
2237 : lb_max => basis_set_b%lmax
2238 : lb_min => basis_set_b%lmin
2239 : npgfb => basis_set_b%npgf
2240 : nsgfb => basis_set_b%nsgf_set
2241 : rpgfb => basis_set_b%pgf_radius
2242 : set_radius_b => basis_set_b%set_radius
2243 : sphi_b => basis_set_b%sphi
2244 : zetb => basis_set_b%zet
2245 : scon_b => basis_set_b%scon
2246 :
2247 : nseta = basis_set_a%nset
2248 : nsetb = basis_set_b%nset
2249 :
2250 : IF (do_symmetric) THEN
2251 : IF (iatom <= jatom) THEN
2252 : irow = iatom
2253 : icol = jatom
2254 : ELSE
2255 : irow = jatom
2256 : icol = iatom
2257 : END IF
2258 : ELSE
2259 : irow = iatom
2260 : icol = jatom
2261 : END IF
2262 :
2263 : NULLIFY (kx_block, ky_block, kz_block)
2264 : CALL dbcsr_get_block_p(matrix_hr(1)%matrix, irow, icol, kx_block, found)
2265 : CPASSERT(found)
2266 : CALL dbcsr_get_block_p(matrix_hr(2)%matrix, irow, icol, ky_block, found)
2267 : CPASSERT(found)
2268 : CALL dbcsr_get_block_p(matrix_hr(3)%matrix, irow, icol, kz_block, found)
2269 : CPASSERT(found)
2270 :
2271 : IF (direction_Or) THEN
2272 : DO ir = 1, 3
2273 : !$OMP CRITICAL(blockadd)
2274 : SELECT CASE (ir)
2275 : CASE (1)
2276 : kx_block(:, :) = kx_block(:, :)*deltaR(ir, jatom)
2277 : CASE (2)
2278 : ky_block(:, :) = ky_block(:, :)*deltaR(ir, jatom)
2279 : CASE (3)
2280 : kz_block(:, :) = kz_block(:, :)*deltaR(ir, jatom)
2281 : END SELECT
2282 : !$OMP END CRITICAL(blockadd)
2283 : END DO
2284 : ELSE
2285 : DO ir = 1, 3
2286 : !$OMP CRITICAL(blockadd)
2287 : SELECT CASE (ir)
2288 : CASE (1)
2289 : kx_block(:, :) = kx_block(:, :)*deltaR(ir, iatom)
2290 : CASE (2)
2291 : ky_block(:, :) = ky_block(:, :)*deltaR(ir, iatom)
2292 : CASE (3)
2293 : kz_block(:, :) = kz_block(:, :)*deltaR(ir, iatom)
2294 : END SELECT
2295 : !$OMP END CRITICAL(blockadd)
2296 : END DO
2297 : END IF
2298 : END DO
2299 : !$OMP END PARALLEL
2300 18 : CALL neighbor_list_iterator_release(nl_iterator)
2301 :
2302 : ! Release work storage
2303 18 : DEALLOCATE (basis_set_list)
2304 :
2305 18 : CALL timestop(handle)
2306 :
2307 36 : END SUBROUTINE hr_mult_by_delta_3d
2308 :
2309 1026 : END MODULE qs_vcd_ao
|
