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