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