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