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