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