(git:691081d)
Loading...
Searching...
No Matches
commutator_rpnl.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!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian
9!> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
10!> \par History
11!> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01]
12!> - full rewrite [jhu, 2009-01-23]
13! **************************************************************************************************
14MODULE commutator_rpnl
15 USE basis_set_types, ONLY: gto_basis_set_p_type,&
16 gto_basis_set_type
17 USE block_p_types, ONLY: block_p_type
18 USE cell_types, ONLY: cell_type
19 USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
20 dbcsr_p_type
21 USE kinds, ONLY: dp
22 USE particle_types, ONLY: particle_type
23 USE qs_kind_types, ONLY: get_qs_kind,&
24 qs_kind_type
25 USE qs_neighbor_list_types, ONLY: get_neighbor_list_set_p,&
26 neighbor_list_set_p_type
27 USE sap_kind_types, ONLY: alist_type,&
28 build_sap_ints,&
29 get_alist,&
30 release_sap_int,&
31 sap_int_type,&
32 sap_sort
33
34!$ USE OMP_LIB, ONLY: omp_lock_kind, &
35!$ omp_init_lock, omp_set_lock, &
36!$ omp_unset_lock, omp_destroy_lock
37
38#include "./base/base_uses.f90"
39
40 IMPLICIT NONE
41
42 PRIVATE
43
44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'commutator_rpnl'
45
46 PUBLIC :: build_com_mom_nl, build_com_nl_mag, build_com_vnl_giao
47
48CONTAINS
49
50! **************************************************************************************************
51!> \brief Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv)
52!> or [rr,Vnl] (matrix_rrv) in AO basis.
53!> Reference point is required for the two latter options
54!> Update: Calculate rxVnlxr (matrix_rvr) and rxrxVnl + Vnlxrxr (matrix_rrv_vrr)
55!> in AO basis. Added in the first place for current correction in
56!> the VG formalism (first order wrt vector potential).
57!> \param qs_kind_set ...
58!> \param sab_all ...
59!> \param sap_ppnl ...
60!> \param eps_ppnl ...
61!> \param particle_set ...
62!> \param cell ...
63!> \param matrix_rv ...
64!> \param matrix_rxrv ...
65!> \param matrix_rrv ...
66!> \param matrix_rvr ...
67!> \param matrix_rrv_vrr ...
68!> \param matrix_r_rxvr ...
69!> \param matrix_rxvr_r ...
70!> \param matrix_r_doublecom ...
71!> \param pseudoatom ...
72!> \param ref_point ...
73! **************************************************************************************************
74 SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, &
75 matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
76
77 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
78 POINTER :: qs_kind_set
79 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
80 INTENT(IN), POINTER :: sab_all, sap_ppnl
81 REAL(kind=dp), INTENT(IN) :: eps_ppnl
82 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
83 POINTER :: particle_set
84 TYPE(cell_type), INTENT(IN), POINTER :: cell
85 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
86 OPTIONAL :: matrix_rv, matrix_rxrv, matrix_rrv, &
87 matrix_rvr, matrix_rrv_vrr
88 TYPE(dbcsr_p_type), DIMENSION(:, :), &
89 INTENT(INOUT), OPTIONAL :: matrix_r_rxvr, matrix_rxvr_r, &
90 matrix_r_doublecom
91 INTEGER, INTENT(in), OPTIONAL :: pseudoatom
92 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
93
94 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_mom_nl'
95 INTEGER, PARAMETER :: i_x = 2, i_xx = 5, i_xy = 6, i_xz = 7, i_y = 3, i_yx = i_xy, i_yy = 8, &
96 i_yz = 9, i_z = 4, i_zx = i_xz, i_zy = i_yz, i_zz = 10
97
98 INTEGER :: handle, i, iab, iac, iatom, ibc, icol, &
99 ikind, ind, ind2, irow, jatom, jkind, &
100 kac, kbc, kkind, na, natom, nb, nkind, &
101 np, order, slot
102 INTEGER, DIMENSION(3) :: cell_b
103 LOGICAL :: asso_r_doublecom, asso_r_rxvr, asso_rrv, asso_rrv_vrr, asso_rv, asso_rvr, &
104 asso_rxrv, asso_rxvr_r, do_symmetric, found, go, my_r_doublecom, my_r_rxvr, my_ref, &
105 my_rrv, my_rrv_vrr, my_rv, my_rvr, my_rxrv, my_rxvr_r, periodic, ppnl_present, trans
106 REAL(kind=dp), DIMENSION(3) :: rab, rf
107 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
108 TYPE(alist_type), POINTER :: alist_ac, alist_bc
109 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_rrv, blocks_rrv_vrr, blocks_rv, &
110 blocks_rvr, blocks_rxrv
111 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_r_doublecom, blocks_r_rxvr, &
112 blocks_rxvr_r
113 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
114 DIMENSION(:) :: basis_set
115 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
116 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
117
118!$ INTEGER(kind=omp_lock_kind), &
119!$ ALLOCATABLE, DIMENSION(:) :: locks
120!$ INTEGER :: lock_num, hash
121!$ INTEGER, PARAMETER :: nlock = 501
122
123 ppnl_present = ASSOCIATED(sap_ppnl)
124 IF (.NOT. ppnl_present) RETURN
125
126 CALL timeset(routinen, handle)
127
128 my_r_doublecom = .false.
129 my_r_rxvr = .false.
130 my_rxvr_r = .false.
131 my_rxrv = .false.
132 my_rrv = .false.
133 my_rv = .false.
134 my_rvr = .false.
135 my_rrv_vrr = .false.
136 IF (PRESENT(matrix_r_doublecom)) my_r_doublecom = .true.
137 IF (PRESENT(matrix_r_rxvr)) my_r_rxvr = .true.
138 IF (PRESENT(matrix_rxvr_r)) my_rxvr_r = .true.
139 IF (PRESENT(matrix_rxrv)) my_rxrv = .true.
140 IF (PRESENT(matrix_rrv)) my_rrv = .true.
141 IF (PRESENT(matrix_rv)) my_rv = .true.
142 IF (PRESENT(matrix_rvr)) my_rvr = .true.
143 IF (PRESENT(matrix_rrv_vrr)) my_rrv_vrr = .true.
144 IF (.NOT. (my_rv .OR. my_rxrv .OR. my_rrv .OR. my_rvr .OR. my_rrv_vrr .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom)) THEN
145 cpabort('No dbcsr matrix provided for commutator calculation!')
146 END IF
147
148 natom = SIZE(particle_set)
149
150 IF (my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) THEN
151 order = 2
152 cpassert(PRESENT(ref_point)) ! need reference point for r x [r,Vnl] and [rr,Vnl]
153 ELSE IF (my_rvr .OR. my_rrv_vrr) THEN
154 order = 2
155 ELSE
156 order = 1
157 END IF
158
159 ! When we want the double commutator [[Vnl, r], r], we also want to fix the pseudoatom
160 IF (my_r_doublecom) THEN
161 cpassert(PRESENT(pseudoatom))
162 END IF
163
164 periodic = any(cell%perd > 0)
165 my_ref = .false.
166 IF (PRESENT(ref_point)) THEN
167 IF (.NOT. periodic) THEN
168 rf = ref_point
169 my_ref = .true.
170 ELSE ! use my_ref = False in periodic case, corresponds to distributed ref point
171 IF (order > 1) THEN
172 cpwarn("Not clear how to define reference point for order > 1 in periodic cells.")
173 END IF
174 END IF
175 END IF
176
177 nkind = SIZE(qs_kind_set)
178
179 !sap_int needs to be shared as multiple threads need to access this
180 NULLIFY (sap_int)
181 ALLOCATE (sap_int(nkind*nkind))
182 DO i = 1, nkind*nkind
183 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
184 sap_int(i)%nalist = 0
185 END DO
186
187 IF (my_ref) THEN
188 ! calculate integrals <a|x^n|p>
189 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=rf, &
190 particle_set=particle_set, cell=cell)
191 ELSE
192 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
193 END IF
194
195 ! *** Set up a sorting index
196 CALL sap_sort(sap_int)
197
198 ALLOCATE (basis_set(nkind))
199 DO ikind = 1, nkind
200 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
201 IF (ASSOCIATED(orb_basis_set)) THEN
202 basis_set(ikind)%gto_basis_set => orb_basis_set
203 ELSE
204 NULLIFY (basis_set(ikind)%gto_basis_set)
205 END IF
206 END DO
207
208 ! *** All integrals needed have been calculated and stored in sap_int
209 ! *** We now calculate the commutator matrix elements
210 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_all, symmetric=do_symmetric)
211
212!$OMP PARALLEL &
213!$OMP DEFAULT (NONE) &
214!$OMP SHARED (basis_set, matrix_rv, matrix_rxrv, matrix_rrv, &
215!$OMP matrix_rvr, matrix_rrv_vrr, matrix_r_doublecom, &
216!$OMP sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
217!$OMP my_rv, my_rxrv, my_rrv, my_rvr, my_rrv_vrr, &
218!$OMP my_r_doublecom, &
219!$OMP matrix_r_rxvr, matrix_rxvr_r, my_r_rxvr, my_rxvr_r, &
220!$OMP pseudoatom, do_symmetric) &
221!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, ind, ind2, &
222!$OMP iab, irow, icol, lock_num, &
223!$OMP blocks_rv, blocks_rxrv, blocks_rrv, blocks_rvr, blocks_rrv_vrr, &
224!$OMP blocks_r_rxvr, blocks_rxvr_r, blocks_r_doublecom, &
225!$OMP found, iac, ibc, alist_ac, alist_bc, &
226!$OMP na, np, nb, kkind, kac, kbc, i, &
227!$OMP go, asso_rv, asso_rxrv, asso_rrv, asso_rvr, asso_rrv_vrr, &
228!$OMP asso_r_rxvr, asso_rxvr_r, asso_r_doublecom, hash, &
229!$OMP acint, achint, bcint, bchint, trans)
230
231!$OMP SINGLE
232!$ ALLOCATE (locks(nlock))
233!$OMP END SINGLE
234
235!$OMP DO
236!$ DO lock_num = 1, nlock
237!$ call omp_init_lock(locks(lock_num))
238!$ END DO
239!$OMP END DO
240
241!$OMP DO SCHEDULE(GUIDED)
242
243 DO slot = 1, sab_all(1)%nl_size
244
245 ikind = sab_all(1)%nlist_task(slot)%ikind
246 jkind = sab_all(1)%nlist_task(slot)%jkind
247 iatom = sab_all(1)%nlist_task(slot)%iatom
248 jatom = sab_all(1)%nlist_task(slot)%jatom
249 cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
250 rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
251
252 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
253 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
254 iab = ikind + nkind*(jkind - 1)
255
256 IF (do_symmetric) THEN
257 IF (iatom <= jatom) THEN
258 irow = iatom
259 icol = jatom
260 ELSE
261 irow = jatom
262 icol = iatom
263 END IF
264 ELSE
265 irow = iatom
266 icol = jatom
267 END IF
268 trans = do_symmetric .AND. (iatom > jatom)
269
270 ! allocate blocks
271 IF (my_rv) THEN
272 ALLOCATE (blocks_rv(3))
273 END IF
274 IF (my_rxrv) THEN
275 ALLOCATE (blocks_rxrv(3))
276 END IF
277 IF (my_rrv) THEN
278 ALLOCATE (blocks_rrv(6))
279 END IF
280 IF (my_rvr) THEN
281 ALLOCATE (blocks_rvr(6))
282 END IF
283 IF (my_rrv_vrr) THEN
284 ALLOCATE (blocks_rrv_vrr(6))
285 END IF
286 IF (my_r_rxvr) THEN
287 ALLOCATE (blocks_r_rxvr(3, 3))
288 END IF
289
290 IF (my_rxvr_r) THEN
291 ALLOCATE (blocks_rxvr_r(3, 3))
292 END IF
293
294 IF (my_r_doublecom) THEN
295 ALLOCATE (blocks_r_doublecom(3, 3))
296 END IF
297
298 ! get blocks
299 IF (my_rv) THEN
300 DO ind = 1, 3
301 CALL dbcsr_get_block_p(matrix_rv(ind)%matrix, irow, icol, blocks_rv(ind)%block, found)
302 END DO
303 END IF
304
305 IF (my_rxrv) THEN
306 DO ind = 1, 3
307 CALL dbcsr_get_block_p(matrix_rxrv(ind)%matrix, irow, icol, blocks_rxrv(ind)%block, found)
308 blocks_rxrv(ind)%block(:, :) = 0._dp
309 END DO
310 END IF
311
312 IF (my_rrv) THEN
313 DO ind = 1, 6
314 CALL dbcsr_get_block_p(matrix_rrv(ind)%matrix, irow, icol, blocks_rrv(ind)%block, found)
315 END DO
316 END IF
317
318 IF (my_rvr) THEN
319 DO ind = 1, 6
320 CALL dbcsr_get_block_p(matrix_rvr(ind)%matrix, irow, icol, blocks_rvr(ind)%block, found)
321 END DO
322 END IF
323
324 IF (my_rrv_vrr) THEN
325 DO ind = 1, 6
326 CALL dbcsr_get_block_p(matrix_rrv_vrr(ind)%matrix, irow, icol, blocks_rrv_vrr(ind)%block, found)
327 END DO
328 END IF
329
330 IF (my_r_rxvr) THEN
331 DO ind = 1, 3
332 DO ind2 = 1, 3
333 CALL dbcsr_get_block_p(matrix_r_rxvr(ind, ind2)%matrix, irow, icol, &
334 blocks_r_rxvr(ind, ind2)%block, found)
335 blocks_r_rxvr(ind, ind2)%block(:, :) = 0._dp
336 END DO
337 END DO
338 END IF
339
340 IF (my_rxvr_r) THEN
341 DO ind = 1, 3
342 DO ind2 = 1, 3
343 CALL dbcsr_get_block_p(matrix_rxvr_r(ind, ind2)%matrix, irow, icol, &
344 blocks_rxvr_r(ind, ind2)%block, found)
345 blocks_rxvr_r(ind, ind2)%block(:, :) = 0._dp
346 END DO
347 END DO
348 END IF
349
350 IF (my_r_doublecom) THEN
351 DO ind = 1, 3
352 DO ind2 = 1, 3
353 CALL dbcsr_get_block_p(matrix_r_doublecom(ind, ind2)%matrix, irow, icol, &
354 blocks_r_doublecom(ind, ind2)%block, found)
355 blocks_r_doublecom(ind, ind2)%block(:, :) = 0._dp
356 END DO
357 END DO
358 END IF
359
360 ! check whether all blocks are associated
361 go = .true.
362 IF (my_rv) THEN
363 asso_rv = (ASSOCIATED(blocks_rv(1)%block) .AND. ASSOCIATED(blocks_rv(2)%block) .AND. &
364 ASSOCIATED(blocks_rv(3)%block))
365 go = go .AND. asso_rv
366 END IF
367
368 IF (my_rxrv) THEN
369 asso_rxrv = (ASSOCIATED(blocks_rxrv(1)%block) .AND. ASSOCIATED(blocks_rxrv(2)%block) .AND. &
370 ASSOCIATED(blocks_rxrv(3)%block))
371 go = go .AND. asso_rxrv
372 END IF
373
374 IF (my_rrv) THEN
375 asso_rrv = (ASSOCIATED(blocks_rrv(1)%block) .AND. ASSOCIATED(blocks_rrv(2)%block) .AND. &
376 ASSOCIATED(blocks_rrv(3)%block) .AND. ASSOCIATED(blocks_rrv(4)%block) .AND. &
377 ASSOCIATED(blocks_rrv(5)%block) .AND. ASSOCIATED(blocks_rrv(6)%block))
378 go = go .AND. asso_rrv
379 END IF
380
381 IF (my_rvr) THEN
382 asso_rvr = (ASSOCIATED(blocks_rvr(1)%block) .AND. ASSOCIATED(blocks_rvr(2)%block) .AND. &
383 ASSOCIATED(blocks_rvr(3)%block) .AND. ASSOCIATED(blocks_rvr(4)%block) .AND. &
384 ASSOCIATED(blocks_rvr(5)%block) .AND. ASSOCIATED(blocks_rvr(6)%block))
385 go = go .AND. asso_rvr
386 END IF
387
388 IF (my_rrv_vrr) THEN
389 asso_rrv_vrr = (ASSOCIATED(blocks_rrv_vrr(1)%block) .AND. ASSOCIATED(blocks_rrv_vrr(2)%block) .AND. &
390 ASSOCIATED(blocks_rrv_vrr(3)%block) .AND. ASSOCIATED(blocks_rrv_vrr(4)%block) .AND. &
391 ASSOCIATED(blocks_rrv_vrr(5)%block) .AND. ASSOCIATED(blocks_rrv_vrr(6)%block))
392 go = go .AND. asso_rrv_vrr
393 END IF
394
395 IF (my_r_rxvr) THEN
396 asso_r_rxvr = .true.
397 DO ind = 1, 3
398 DO ind2 = 1, 3
399 asso_r_rxvr = asso_r_rxvr .AND. ASSOCIATED(blocks_r_rxvr(ind, ind2)%block)
400 END DO
401 END DO
402 go = go .AND. asso_r_rxvr
403 END IF
404
405 IF (my_rxvr_r) THEN
406 asso_rxvr_r = .true.
407 DO ind = 1, 3
408 DO ind2 = 1, 3
409 asso_rxvr_r = asso_rxvr_r .AND. ASSOCIATED(blocks_rxvr_r(ind, ind2)%block)
410 END DO
411 END DO
412 go = go .AND. asso_rxvr_r
413 END IF
414
415 IF (my_r_doublecom) THEN
416 asso_r_doublecom = .true.
417 DO ind = 1, 3
418 DO ind2 = 1, 3
419 asso_r_doublecom = asso_r_doublecom .AND. ASSOCIATED(blocks_r_doublecom(ind, ind2)%block)
420 END DO
421 END DO
422 go = go .AND. asso_r_doublecom
423 END IF
424
425 ! loop over all kinds for projector atom
426 ! < iatom | katom > h < katom | jatom >
427 IF (go) THEN
428 DO kkind = 1, nkind
429 iac = ikind + nkind*(kkind - 1)
430 ibc = jkind + nkind*(kkind - 1)
431 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
432 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
433 CALL get_alist(sap_int(iac), alist_ac, iatom)
434 CALL get_alist(sap_int(ibc), alist_bc, jatom)
435 IF (.NOT. ASSOCIATED(alist_ac)) cycle
436 IF (.NOT. ASSOCIATED(alist_bc)) cycle
437 DO kac = 1, alist_ac%nclist
438 DO kbc = 1, alist_bc%nclist
439 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
440 IF (PRESENT(pseudoatom)) THEN
441 IF (alist_ac%clist(kac)%catom /= pseudoatom) cycle
442 END IF
443
444 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
445 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
446 acint => alist_ac%clist(kac)%acint
447 bcint => alist_bc%clist(kbc)%acint
448 achint => alist_ac%clist(kac)%achint
449 bchint => alist_bc%clist(kbc)%achint
450 na = SIZE(acint, 1)
451 np = SIZE(acint, 2)
452 nb = SIZE(bcint, 1)
453!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
454!$ CALL omp_set_lock(locks(hash))
455 IF (my_rv) THEN
456 ! r*Vnl
457 ! with LAPACK
458 ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 2), na, &
459 ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(1)%block, SIZE(blocks_rv(1)%block, 1)) ! xV
460 ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 3), na, &
461 ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(2)%block, SIZE(blocks_rv(2)%block, 1)) ! yV
462 ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 4), na, &
463 ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(3)%block, SIZE(blocks_rv(3)%block, 1)) ! zV
464 IF (.NOT. trans) THEN
465 ! with MATMUL
466 blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
467 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1))) ! xV
468 blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) + &
469 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1))) ! yV
470 blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) + &
471 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1))) ! zV
472 ELSE
473 blocks_rv(1)%block(1:nb, 1:na) = blocks_rv(1)%block(1:nb, 1:na) + &
474 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))
475 blocks_rv(2)%block(1:nb, 1:na) = blocks_rv(2)%block(1:nb, 1:na) + &
476 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))
477 blocks_rv(3)%block(1:nb, 1:na) = blocks_rv(3)%block(1:nb, 1:na) + &
478 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))
479 END IF
480 ! -Vnl r
481 ! with LAPACK
482 ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
483 ! bcint(1, 1, 2), nb, 1.0_dp, blocks_rv(1)%block, SIZE(blocks_rv(1)%block, 1)) ! -Vx
484 ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
485 ! bcint(1, 1, 3), nb, 1.0_dp, blocks_rv(2)%block, SIZE(blocks_rv(2)%block, 1)) ! -Vy
486 ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
487 ! bcint(1, 1, 4), nb, 1.0_dp, blocks_rv(3)%block, SIZE(blocks_rv(3)%block, 1)) ! -Vz
488 ! with MATMUL
489 IF (.NOT. trans) THEN
490 blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) - &
491 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) ! -Vx
492 blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) - &
493 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) ! -Vy
494 blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) - &
495 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) ! -Vz
496 ELSE
497 blocks_rv(1)%block(1:nb, 1:na) = blocks_rv(1)%block(1:nb, 1:na) - &
498 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2)))
499 blocks_rv(2)%block(1:nb, 1:na) = blocks_rv(2)%block(1:nb, 1:na) - &
500 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3)))
501 blocks_rv(3)%block(1:nb, 1:na) = blocks_rv(3)%block(1:nb, 1:na) - &
502 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4)))
503 END IF
504
505 END IF
506
507 IF (my_rxrv) THEN
508 ! x-component (y [z,Vnl] - z [y, Vnl])
509 IF (iatom <= jatom) THEN
510 ! yzV
511 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
512 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
513 ! -yVz
514 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
515 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4)))
516 ! -zyV
517 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
518 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
519 ! zVy
520 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
521 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 3)))
522 ELSE
523 ! yzV
524 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
525 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
526 ! -yVz
527 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
528 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4)))
529 ! -zyV
530 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
531 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
532 ! zVy
533 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
534 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 3)))
535 END IF
536
537 ! y-component (z [x,Vnl] - x [z, Vnl])
538 IF (iatom <= jatom) THEN
539 ! zxV
540 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
541 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
542 ! -zVx
543 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
544 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 2)))
545 ! -xzV
546 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
547 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
548 ! xVz
549 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
550 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4)))
551 ELSE
552 ! zxV
553 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
554 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
555 ! -zVx
556 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
557 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 2)))
558 ! -xzV
559 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
560 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
561 ! xVz
562 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
563 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4)))
564 END IF
565
566 ! z-component (x [y,Vnl] - y [x, Vnl])
567 IF (iatom <= jatom) THEN
568 ! xyV
569 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
570 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
571 ! -xVy
572 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
573 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3)))
574 ! -yxV
575 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
576 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
577 ! zVx
578 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
579 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 2)))
580 ELSE
581 ! xyV
582 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
583 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
584 ! -xVy
585 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
586 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3)))
587 ! -yxV
588 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
589 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
590 ! zVx
591 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
592 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 2)))
593 END IF
594 END IF
595
596 IF (my_rrv) THEN
597 ! r_alpha * r_beta * Vnl
598 IF (iatom <= jatom) THEN
599 ! xxV
600 blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) + &
601 matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1)))
602 ! xyV
603 blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) + &
604 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
605 ! xzV
606 blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) + &
607 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
608 ! yyV
609 blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) + &
610 matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1)))
611 ! yzV
612 blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) + &
613 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
614 ! zzV
615 blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) + &
616 matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1)))
617 ELSE
618 ! xxV
619 blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) + &
620 matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1)))
621 ! xyV
622 blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) + &
623 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
624 ! xzV
625 blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) + &
626 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
627 ! yyV
628 blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) + &
629 matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1)))
630 ! yzV
631 blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) + &
632 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
633 ! zzV
634 blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) + &
635 matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1)))
636 END IF
637
638 ! - Vnl * r_alpha * r_beta
639 IF (iatom <= jatom) THEN
640 ! -Vxx
641 blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) - &
642 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5)))
643 ! -Vxy
644 blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) - &
645 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6)))
646 ! -Vxz
647 blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) - &
648 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7)))
649 ! -Vyy
650 blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) - &
651 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8)))
652 ! -Vyz
653 blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) - &
654 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9)))
655 ! -Vzz
656 blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) - &
657 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10)))
658 ELSE
659 ! -Vxx
660 blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) - &
661 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5)))
662 ! -Vxy
663 blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) - &
664 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6)))
665 ! -Vxz
666 blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) - &
667 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7)))
668 ! -Vyy
669 blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) - &
670 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8)))
671 ! -Vyz
672 blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) - &
673 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9)))
674 ! -Vzz
675 blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) - &
676 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10)))
677 END IF
678 END IF
679
680 IF (my_rvr) THEN
681 ! r_alpha * Vnl * r_beta
682 IF (iatom <= jatom) THEN
683 ! xVx
684 blocks_rvr(1)%block(1:na, 1:nb) = blocks_rvr(1)%block(1:na, 1:nb) + &
685 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 2)))
686 ! xVy
687 blocks_rvr(2)%block(1:na, 1:nb) = blocks_rvr(2)%block(1:na, 1:nb) + &
688 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3)))
689 ! xVz
690 blocks_rvr(3)%block(1:na, 1:nb) = blocks_rvr(3)%block(1:na, 1:nb) + &
691 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4)))
692 ! yVy
693 blocks_rvr(4)%block(1:na, 1:nb) = blocks_rvr(4)%block(1:na, 1:nb) + &
694 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 3)))
695 ! yVz
696 blocks_rvr(5)%block(1:na, 1:nb) = blocks_rvr(5)%block(1:na, 1:nb) + &
697 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4)))
698 ! zVz
699 blocks_rvr(6)%block(1:na, 1:nb) = blocks_rvr(6)%block(1:na, 1:nb) + &
700 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 4)))
701 ELSE
702 ! xVx
703 blocks_rvr(1)%block(1:nb, 1:na) = blocks_rvr(1)%block(1:nb, 1:na) + &
704 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 2)))
705 ! xVy
706 blocks_rvr(2)%block(1:nb, 1:na) = blocks_rvr(2)%block(1:nb, 1:na) + &
707 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3)))
708 ! xVz
709 blocks_rvr(3)%block(1:nb, 1:na) = blocks_rvr(3)%block(1:nb, 1:na) + &
710 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4)))
711 ! yVy
712 blocks_rvr(4)%block(1:nb, 1:na) = blocks_rvr(4)%block(1:nb, 1:na) + &
713 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 3)))
714 ! yVz
715 blocks_rvr(5)%block(1:nb, 1:na) = blocks_rvr(5)%block(1:nb, 1:na) + &
716 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4)))
717 ! zVz
718 blocks_rvr(6)%block(1:nb, 1:na) = blocks_rvr(6)%block(1:nb, 1:na) + &
719 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 4)))
720 END IF
721 END IF
722
723 IF (my_rrv_vrr) THEN
724 ! r_alpha * r_beta * Vnl
725 IF (iatom <= jatom) THEN
726 ! xxV
727 blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
728 matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1)))
729 ! xyV
730 blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
731 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
732 ! xzV
733 blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
734 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
735 ! yyV
736 blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
737 matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1)))
738 ! yzV
739 blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
740 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
741 ! zzV
742 blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
743 matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1)))
744 ELSE
745 ! xxV
746 blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
747 matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1)))
748 ! xyV
749 blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
750 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
751 ! xzV
752 blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
753 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
754 ! yyV
755 blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
756 matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1)))
757 ! yzV
758 blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
759 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
760 ! zzV
761 blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
762 matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1)))
763 END IF
764 ! + Vnl * r_alpha * r_beta
765 IF (iatom <= jatom) THEN
766 ! +Vxx
767 blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
768 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5)))
769 ! +Vxy
770 blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
771 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6)))
772 ! +Vxz
773 blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
774 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7)))
775 ! +Vyy
776 blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
777 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8)))
778 ! +Vyz
779 blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
780 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9)))
781 ! +Vzz
782 blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
783 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10)))
784 ELSE
785 ! +Vxx
786 blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
787 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5)))
788 ! +Vxy
789 blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
790 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6)))
791 ! +Vxz
792 blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
793 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7)))
794 ! +Vyy
795 blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
796 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8)))
797 ! +Vyz
798 blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
799 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9)))
800 ! +Vzz
801 blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
802 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10)))
803 END IF
804 END IF
805
806 ! The indices are stored in i_1, i_x, ..., i_zzz
807
808 ! TODO: is this set to zero before?
809 IF (my_r_rxvr) THEN
810 ! beta = 1
811 ! matrix_r_rxvr(x, x) = x * y * V_nl * z - x * z * V_nl * y
812 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
813 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) + &
814 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
815 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
816 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) - &
817 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
818
819 ! matrix_r_rxvr(y, x) = x * z * V_nl * x - x * x * V_nl * z
820 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
821 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) + &
822 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_x)))
823 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
824 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) - &
825 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
826
827 ! matrix_r_rxvr(z, x) = x * x * V_nl * y - x * y * V_nl * x
828 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
829 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) + &
830 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
831 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
832 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) - &
833 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_x)))
834
835 ! beta = 2
836 ! matrix_r_rxvr(x, y) = y * y * V_nl * z - y * z * V_nl * y
837 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
838 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) + &
839 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
840 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
841 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) - &
842 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_y)))
843
844 ! matrix_r_rxvr(y, y) = y * z * V_nl * x - y * x * V_nl * z
845 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
846 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) + &
847 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
848 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
849 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) - &
850 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
851
852 ! matrix_r_rxvr(z, y) = y * x * V_nl * y - y * y * V_nl * x
853 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
854 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) + &
855 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_y)))
856 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
857 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) - &
858 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
859
860 ! beta = 3
861 ! matrix_r_rxvr(x, z) = z * y * V_nl * z - z * z * V_nl * y
862 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
863 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) + &
864 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_z)))
865 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
866 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) - &
867 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
868
869 ! matrix_r_rxvr(y, z) = z * z * V_nl * x - z * x * V_nl * z
870 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
871 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) + &
872 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
873 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
874 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) - &
875 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_z)))
876
877 ! matrix_r_rxvr(z, z) = z * x * V_nl * y - z * y * V_nl * x
878 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
879 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) + &
880 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
881 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
882 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) - &
883 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
884
885 END IF ! my_r_rxvr
886
887 ! The indices are stored in i_1, i_x, ..., i_zzz
888 ! This will put into blocks_rxvr_r
889 ! matrix_rxvr_r(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
890 ! r_gamma * V_nl * r_delta * r_beta
891 IF (my_rxvr_r) THEN
892 ! beta = 1
893 ! matrix_rxvr_r(x, x) = yV zx - zV yx
894 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
895 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) + &
896 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
897 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
898 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) - &
899 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
900
901 ! matrix_rxvr_r(y, x) = zV xx - xV zx
902 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
903 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) + &
904 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
905 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
906 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) - &
907 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zx)))
908
909 ! matrix_rxvr_r(z, x) = xV yx - yV xx
910 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
911 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) + &
912 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yx)))
913 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
914 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) - &
915 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
916
917 ! beta = 2
918 ! matrix_rxvr_r(x, y) = yV zy - zV yy
919 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
920 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) + &
921 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zy)))
922 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
923 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) - &
924 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
925
926 ! matrix_rxvr_r(y, y) = zV xy - xV zy
927 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
928 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) + &
929 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
930 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
931 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) - &
932 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
933
934 ! matrix_rxvr_r(z, y) = xV yy - yV xy
935 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
936 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) + &
937 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
938 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
939 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) - &
940 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xy)))
941
942 ! beta = 3
943 ! matrix_rxvr_r(x, z) = yV zz - zV yz
944 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
945 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) + &
946 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
947 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
948 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) - &
949 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yz)))
950
951 ! matrix_rxvr_r(y, z) = zV xz - xV zz
952 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
953 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) + &
954 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xz)))
955 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
956 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) - &
957 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
958
959 ! matrix_rxvr_r(z, z) = xV yz - yV xz
960 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
961 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) + &
962 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
963 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
964 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) - &
965 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
966
967 END IF ! my_rxvr_r
968
969 ! matrix_r_doublecom(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
970 ! gamma V^pseudoatom beta delta - gamma beta V^pseudoatom delta
971
972 IF (my_r_doublecom) THEN
973 ! beta = 1
974 ! matrix_r_doublecom(x, x) = yV xz - zV xy - yxV z + zxV y
975 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
976 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
977 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
978 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
979 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
980 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
981 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
982 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
983 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
984 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
985 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
986 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
987
988 ! matrix_r_doublecom(y, x) = zV xx - xV xz - zxV x + xxV z
989 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
990 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
991 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
992 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
993 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
994 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xz)))
995 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
996 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
997 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_x)))
998 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
999 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
1000 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
1001
1002 ! matrix_r_doublecom(z, x) = xV xy - yV xx - xxV y + yxV x
1003 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1004 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1005 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xy)))
1006 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1007 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1008 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
1009 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1010 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1011 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
1012 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1013 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1014 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_x)))
1015
1016 ! beta = 2
1017 ! matrix_r_doublecom(x, y) = yV yz - zV yy - yyV z + zyV y
1018 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1019 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1020 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yz)))
1021 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1022 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1023 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
1024 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1025 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1026 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
1027 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1028 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1029 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_y)))
1030
1031 ! matrix_r_doublecom(y, y) = zV yx - xV yz - zyV x + xyV z
1032 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1033 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1034 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
1035 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1036 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1037 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
1038 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1039 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1040 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
1041 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1042 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1043 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
1044
1045 ! matrix_r_doublecom(z, y) = xV yy - yV yx - xyV y + yyV x
1046 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1047 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1048 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
1049 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1050 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1051 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yx)))
1052 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1053 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1054 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_y)))
1055 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1056 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1057 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
1058
1059 ! beta = 3
1060 ! matrix_r_doublecom(x, z) = yV zz - zV zy - yzV z + zzV y
1061 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1062 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1063 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
1064 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1065 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1066 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zy)))
1067 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1068 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1069 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_z)))
1070 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1071 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1072 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
1073
1074 ! matrix_r_doublecom(y, z) = zV zx - xV zz - zzV x + xzV z
1075 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1076 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1077 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zx)))
1078 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1079 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1080 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
1081 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1082 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1083 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
1084 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1085 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1086 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_z)))
1087
1088 ! matrix_r_doublecom(z, z) = xV zy - yV zx - xzV y + yzV x
1089 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1090 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1091 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
1092 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1093 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1094 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
1095 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1096 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1097 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
1098 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1099 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1100 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
1101
1102 END IF ! my_r_doublecom
1103!$ CALL omp_unset_lock(locks(hash))
1104 EXIT ! We have found a match and there can be only one single match
1105 END IF
1106 END DO
1107 END DO
1108 END DO
1109 END IF
1110 IF (my_rv) THEN
1111 DO ind = 1, 3
1112 NULLIFY (blocks_rv(ind)%block)
1113 END DO
1114 DEALLOCATE (blocks_rv)
1115 END IF
1116 IF (my_rxrv) THEN
1117 DO ind = 1, 3
1118 NULLIFY (blocks_rxrv(ind)%block)
1119 END DO
1120 DEALLOCATE (blocks_rxrv)
1121 END IF
1122 IF (my_rrv) THEN
1123 DO ind = 1, 6
1124 NULLIFY (blocks_rrv(ind)%block)
1125 END DO
1126 DEALLOCATE (blocks_rrv)
1127 END IF
1128 IF (my_rvr) THEN
1129 DO ind = 1, 6
1130 NULLIFY (blocks_rvr(ind)%block)
1131 END DO
1132 DEALLOCATE (blocks_rvr)
1133 END IF
1134 IF (my_rrv_vrr) THEN
1135 DO ind = 1, 6
1136 NULLIFY (blocks_rrv_vrr(ind)%block)
1137 END DO
1138 DEALLOCATE (blocks_rrv_vrr)
1139 END IF
1140 IF (my_r_rxvr) THEN
1141 DO ind = 1, 3
1142 DO ind2 = 1, 3
1143 NULLIFY (blocks_r_rxvr(ind, ind2)%block)
1144 END DO
1145 END DO
1146 DEALLOCATE (blocks_r_rxvr)
1147 END IF
1148 IF (my_rxvr_r) THEN
1149 DO ind = 1, 3
1150 DO ind2 = 1, 3
1151 NULLIFY (blocks_rxvr_r(ind, ind2)%block)
1152 END DO
1153 END DO
1154 DEALLOCATE (blocks_rxvr_r)
1155 END IF
1156 IF (my_r_doublecom) THEN
1157 DO ind = 1, 3
1158 DO ind2 = 1, 3
1159 NULLIFY (blocks_r_doublecom(ind, ind2)%block)
1160 END DO
1161 END DO
1162 DEALLOCATE (blocks_r_doublecom)
1163 END IF
1164 END DO
1165
1166!$OMP DO
1167!$ DO lock_num = 1, nlock
1168!$ call omp_destroy_lock(locks(lock_num))
1169!$ END DO
1170!$OMP END DO
1171
1172!$OMP SINGLE
1173!$ DEALLOCATE (locks)
1174!$OMP END SINGLE NOWAIT
1175
1176!$OMP END PARALLEL
1177
1178 CALL release_sap_int(sap_int)
1179
1180 DEALLOCATE (basis_set)
1181
1182 CALL timestop(handle)
1183
1184 END SUBROUTINE build_com_mom_nl
1185
1186! **************************************************************************************************
1187!> \brief calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
1188!> \param qs_kind_set ...
1189!> \param sab_all ...
1190!> \param sap_ppnl ...
1191!> \param eps_ppnl ...
1192!> \param particle_set ...
1193!> \param matrix_mag_nl ...
1194!> \param refpoint ...
1195!> \param cell ...
1196! **************************************************************************************************
1197 SUBROUTINE build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
1198
1199 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1200 POINTER :: qs_kind_set
1201 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1202 INTENT(IN), POINTER :: sab_all, sap_ppnl
1203 REAL(kind=dp), INTENT(IN) :: eps_ppnl
1204 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1205 POINTER :: particle_set
1206 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1207 POINTER :: matrix_mag_nl
1208 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
1209 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1210
1211 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_nl_mag'
1212
1213 INTEGER :: handle, iab, iac, iatom, ibc, icol, &
1214 ikind, ind, irow, jatom, jkind, kac, &
1215 kbc, kkind, na, natom, nb, nkind, np, &
1216 order, slot
1217 INTEGER, DIMENSION(3) :: cell_b
1218 LOGICAL :: found, go, my_ref, ppnl_present
1219 REAL(kind=dp), DIMENSION(3) :: r_b, r_ps, rab
1220 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1221 TYPE(alist_type), POINTER :: alist_ac, alist_bc
1222 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1223 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1224 DIMENSION(:) :: basis_set
1225 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1226 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1227
1228!$ INTEGER(kind=omp_lock_kind), &
1229!$ ALLOCATABLE, DIMENSION(:) :: locks
1230!$ INTEGER :: lock_num, hash
1231!$ INTEGER, PARAMETER :: nlock = 501
1232
1233 ppnl_present = ASSOCIATED(sap_ppnl)
1234 IF (.NOT. ppnl_present) RETURN
1235
1236 CALL timeset(routinen, handle)
1237
1238 my_ref = .false.
1239 IF (PRESENT(refpoint)) THEN
1240 my_ref = .true.
1241 cpassert(PRESENT(cell))
1242 END IF
1243
1244 natom = SIZE(particle_set)
1245 nkind = SIZE(qs_kind_set)
1246
1247 ! allocate integral storage
1248 NULLIFY (sap_int)
1249 ALLOCATE (sap_int(nkind*nkind))
1250 DO ind = 1, nkind*nkind
1251 NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1252 sap_int(ind)%nalist = 0
1253 END DO
1254
1255 ! build integrals over GTO + projector functions, refpoint actually
1256 order = 1 ! only need first moments (x, y, z)
1257 ! refpoint actually does not matter in this case, i. e. (order = 1 .and. commutator)
1258 IF (my_ref) THEN
1259 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=refpoint, &
1260 particle_set=particle_set, cell=cell)
1261 ELSE
1262 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1263 END IF
1264
1265 CALL sap_sort(sap_int)
1266
1267 ! get access to basis sets
1268 ALLOCATE (basis_set(nkind))
1269 DO ikind = 1, nkind
1270 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1271 IF (ASSOCIATED(orb_basis_set)) THEN
1272 basis_set(ikind)%gto_basis_set => orb_basis_set
1273 ELSE
1274 NULLIFY (basis_set(ikind)%gto_basis_set)
1275 END IF
1276 END DO
1277
1278!$OMP PARALLEL &
1279!$OMP DEFAULT (NONE) &
1280!$OMP SHARED (basis_set, matrix_mag_nl, sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
1281!$OMP particle_set, my_ref, refpoint) &
1282!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, lock_num, &
1283!$OMP iab, irow, icol, blocks_mag, r_ps, r_b, go, hash, &
1284!$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
1285!$OMP achint, bchint, na, np, nb, kkind, kac, kbc)
1286
1287!$OMP SINGLE
1288!$ ALLOCATE (locks(nlock))
1289!$OMP END SINGLE
1290
1291!$OMP DO
1292!$ DO lock_num = 1, nlock
1293!$ call omp_init_lock(locks(lock_num))
1294!$ END DO
1295!$OMP END DO
1296
1297!$OMP DO SCHEDULE(GUIDED)
1298 DO slot = 1, sab_all(1)%nl_size
1299 ! get indices
1300 ikind = sab_all(1)%nlist_task(slot)%ikind
1301 jkind = sab_all(1)%nlist_task(slot)%jkind
1302 iatom = sab_all(1)%nlist_task(slot)%iatom
1303 jatom = sab_all(1)%nlist_task(slot)%jatom
1304 cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1305 rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1306
1307 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1308 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1309 iab = ikind + nkind*(jkind - 1)
1310
1311 IF (iatom <= jatom) THEN
1312 irow = iatom
1313 icol = jatom
1314 ELSE
1315 irow = jatom
1316 icol = iatom
1317 END IF
1318
1319 ! get blocks
1320 ALLOCATE (blocks_mag(3))
1321 DO ind = 1, 3
1322 CALL dbcsr_get_block_p(matrix_mag_nl(ind)%matrix, irow, icol, blocks_mag(ind)%block, found)
1323 END DO
1324
1325 go = (ASSOCIATED(blocks_mag(1)%block) .AND. ASSOCIATED(blocks_mag(2)%block) .AND. ASSOCIATED(blocks_mag(3)%block))
1326
1327 IF (go) THEN
1328 DO kkind = 1, nkind
1329 iac = ikind + nkind*(kkind - 1)
1330 ibc = jkind + nkind*(kkind - 1)
1331 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1332 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1333 CALL get_alist(sap_int(iac), alist_ac, iatom)
1334 CALL get_alist(sap_int(ibc), alist_bc, jatom)
1335 IF (.NOT. ASSOCIATED(alist_ac)) cycle
1336 IF (.NOT. ASSOCIATED(alist_bc)) cycle
1337 DO kac = 1, alist_ac%nclist
1338 DO kbc = 1, alist_bc%nclist
1339 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1340 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1341 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1342
1343 acint => alist_ac%clist(kac)%acint
1344 bcint => alist_bc%clist(kbc)%acint
1345 achint => alist_ac%clist(kac)%achint
1346 bchint => alist_bc%clist(kbc)%achint
1347 na = SIZE(acint, 1)
1348 np = SIZE(acint, 2)
1349 nb = SIZE(bcint, 1)
1350 ! Position of the pseudized atom
1351 r_ps = particle_set(alist_ac%clist(kac)%catom)%r
1352 r_b = refpoint
1353
1354!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1355!$ CALL omp_set_lock(locks(hash))
1356 ! assemble integrals
1357 IF (iatom <= jatom) THEN
1358 blocks_mag(1)%block(1:na, 1:nb) = blocks_mag(1)%block(1:na, 1:nb) + &
1359 (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1360 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) & ! R_y [V_nl, z]
1361 - (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1362 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) ! - R_z [V_nl, y]
1363 blocks_mag(2)%block(1:na, 1:nb) = blocks_mag(2)%block(1:na, 1:nb) + &
1364 (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1365 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) & ! R_z [V_nl, x]
1366 - (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1367 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) ! - R_x [V_nl, z]
1368 blocks_mag(3)%block(1:na, 1:nb) = blocks_mag(3)%block(1:na, 1:nb) + &
1369 (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1370 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) & ! R_x [V_nl, y]
1371 - (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1372 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) ! - R_y [V_nl, x]
1373 ELSE
1374 blocks_mag(1)%block(1:nb, 1:na) = blocks_mag(1)%block(1:nb, 1:na) + &
1375 (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1376 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) & ! R_y [V_nl, z]
1377 - (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1378 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) ! - R_z [V_nl, y]
1379 blocks_mag(2)%block(1:nb, 1:na) = blocks_mag(2)%block(1:nb, 1:na) + &
1380 (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1381 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) & ! R_z [V_nl, x]
1382 - (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1383 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) ! - R_x [V_nl, z]
1384 blocks_mag(3)%block(1:nb, 1:na) = blocks_mag(3)%block(1:nb, 1:na) + &
1385 (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1386 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) & ! R_x [V_nl, y]
1387 - (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1388 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) ! - R_y [V_nl, x]
1389 END IF
1390!$ CALL omp_unset_lock(locks(hash))
1391 EXIT ! We have found a match and there can be only one single match
1392 END IF
1393 END DO
1394 END DO
1395 END DO
1396 END IF
1397
1398 DO ind = 1, 3
1399 NULLIFY (blocks_mag(ind)%block)
1400 END DO
1401 DEALLOCATE (blocks_mag)
1402 END DO
1403
1404!$OMP DO
1405!$ DO lock_num = 1, nlock
1406!$ call omp_destroy_lock(locks(lock_num))
1407!$ END DO
1408!$OMP END DO
1409
1410!$OMP SINGLE
1411!$ DEALLOCATE (locks)
1412!$OMP END SINGLE NOWAIT
1413
1414!$OMP END PARALLEL
1415
1416 DEALLOCATE (basis_set)
1417 CALL release_sap_int(sap_int)
1418
1419 CALL timestop(handle)
1420
1421 END SUBROUTINE build_com_nl_mag
1422
1423! **************************************************************************************************
1424!> \brief Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs
1425!> \param qs_kind_set ...
1426!> \param sab_all ...
1427!> \param sap_ppnl ...
1428!> \param eps_ppnl ...
1429!> \param particle_set ...
1430!> \param matrix_rv ...
1431!> \param ref_point ...
1432!> \param cell ...
1433!> \param direction_Or If set to true: calculate Vnl * r_delta
1434!> Otherwise calculate r_delta * Vnl
1435! **************************************************************************************************
1436 SUBROUTINE build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, &
1437 matrix_rv, ref_point, cell, direction_Or)
1438
1439 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1440 POINTER :: qs_kind_set
1441 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1442 INTENT(IN), POINTER :: sab_all, sap_ppnl
1443 REAL(kind=dp), INTENT(IN) :: eps_ppnl
1444 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1445 POINTER :: particle_set
1446 TYPE(dbcsr_p_type), DIMENSION(:, :), &
1447 INTENT(INOUT), OPTIONAL, POINTER :: matrix_rv
1448 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
1449 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1450 LOGICAL :: direction_or
1451
1452 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_vnl_giao'
1453 INTEGER, PARAMETER :: i_1 = 1
1454
1455 INTEGER :: delta, gamma, handle, i, iab, iac, &
1456 iatom, ibc, icol, ikind, irow, j, &
1457 jatom, jkind, kac, kbc, kkind, na, &
1458 natom, nb, nkind, np, order, slot
1459 INTEGER, DIMENSION(3) :: cell_b
1460 LOGICAL :: found, my_ref, ppnl_present
1461 REAL(kind=dp), DIMENSION(3) :: rab, rf
1462 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1463 TYPE(alist_type), POINTER :: alist_ac, alist_bc
1464 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1465 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1466 DIMENSION(:) :: basis_set
1467 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1468 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1469
1470!$ INTEGER(kind=omp_lock_kind), &
1471!$ ALLOCATABLE, DIMENSION(:) :: locks
1472!$ INTEGER :: lock_num, hash
1473!$ INTEGER, PARAMETER :: nlock = 501
1474
1475 ppnl_present = ASSOCIATED(sap_ppnl)
1476 IF (.NOT. ppnl_present) RETURN
1477
1478 CALL timeset(routinen, handle)
1479
1480 natom = SIZE(particle_set)
1481
1482 my_ref = .false.
1483 IF (PRESENT(ref_point)) THEN
1484 cpassert(PRESENT(cell)) ! need cell as well if refpoint is provided
1485 rf = ref_point
1486 my_ref = .true.
1487 END IF
1488
1489 nkind = SIZE(qs_kind_set)
1490
1491 ! sap_int needs to be shared as multiple threads need to access this
1492 NULLIFY (sap_int)
1493 ALLOCATE (sap_int(nkind*nkind))
1494 DO i = 1, nkind*nkind
1495 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1496 sap_int(i)%nalist = 0
1497 END DO
1498
1499 order = 1
1500 IF (my_ref) THEN
1501 ! calculate integrals <a|x^n|p>
1502 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=rf, &
1503 particle_set=particle_set, cell=cell)
1504 ELSE
1505 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1506 END IF
1507
1508 ! *** Set up a sorting index
1509 CALL sap_sort(sap_int)
1510
1511 ALLOCATE (basis_set(nkind))
1512 DO ikind = 1, nkind
1513 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1514 IF (ASSOCIATED(orb_basis_set)) THEN
1515 basis_set(ikind)%gto_basis_set => orb_basis_set
1516 ELSE
1517 NULLIFY (basis_set(ikind)%gto_basis_set)
1518 END IF
1519 END DO
1520
1521 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_all)
1522 ! *** All integrals needed have been calculated and stored in sap_int
1523 ! *** We now calculate the commutator matrix elements
1524
1525!$OMP PARALLEL &
1526!$OMP DEFAULT (NONE) &
1527!$OMP SHARED (basis_set, matrix_rv, &
1528!$OMP sap_int, nkind, eps_ppnl, locks, sab_all, &
1529!$OMP particle_set, direction_Or) &
1530!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1531!$OMP iab, irow, icol, blocks_rv, &
1532!$OMP found, iac, ibc, alist_ac, alist_bc, &
1533!$OMP na, np, nb, kkind, kac, kbc, i, lock_num, &
1534!$OMP hash, natom, delta, gamma, achint, bchint, acint, bcint)
1535
1536!$OMP SINGLE
1537!$ ALLOCATE (locks(nlock))
1538!$OMP END SINGLE
1539
1540!$OMP DO
1541!$ DO lock_num = 1, nlock
1542!$ call omp_init_lock(locks(lock_num))
1543!$ END DO
1544!$OMP END DO
1545
1546!$OMP DO SCHEDULE(GUIDED)
1547
1548 DO slot = 1, sab_all(1)%nl_size
1549
1550 ikind = sab_all(1)%nlist_task(slot)%ikind
1551 jkind = sab_all(1)%nlist_task(slot)%jkind
1552 iatom = sab_all(1)%nlist_task(slot)%iatom
1553 jatom = sab_all(1)%nlist_task(slot)%jatom
1554 cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1555 rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1556
1557 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1558 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1559 iab = ikind + nkind*(jkind - 1)
1560
1561 irow = iatom
1562 icol = jatom
1563
1564 ! allocate blocks
1565 ALLOCATE (blocks_rv(3, 3))
1566
1567 ! get blocks
1568 DO i = 1, 3
1569 DO j = 1, 3
1570 CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, &
1571 blocks_rv(i, j)%block, found)
1572 blocks_rv(i, j)%block(:, :) = 0.0_dp
1573 cpassert(found)
1574 END DO
1575 END DO
1576
1577 ! loop over all kinds for projector atom
1578 ! < iatom | katom > h < katom | jatom >
1579 DO kkind = 1, nkind
1580 iac = ikind + nkind*(kkind - 1)
1581 ibc = jkind + nkind*(kkind - 1)
1582 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1583 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1584 CALL get_alist(sap_int(iac), alist_ac, iatom)
1585 CALL get_alist(sap_int(ibc), alist_bc, jatom)
1586 IF (.NOT. ASSOCIATED(alist_ac)) cycle
1587 IF (.NOT. ASSOCIATED(alist_bc)) cycle
1588 DO kac = 1, alist_ac%nclist
1589 DO kbc = 1, alist_bc%nclist
1590 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1591
1592 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1593 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1594 acint => alist_ac%clist(kac)%acint
1595 bcint => alist_bc%clist(kbc)%acint
1596 achint => alist_ac%clist(kac)%achint
1597 bchint => alist_bc%clist(kbc)%achint
1598 na = SIZE(acint, 1)
1599 np = SIZE(acint, 2)
1600 nb = SIZE(bcint, 1)
1601!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1602!$ CALL omp_set_lock(locks(hash))
1603
1604 !! The atom index is alist_ac%clist(kac)%catom
1605 ! The coordinate is particle_set(alist_ac%clist(kac)%catom)%r(:)
1606 IF (direction_or) THEN ! V * r_delta * (R^eta_gamma - R^nu_gamma)
1607 DO delta = 1, 3
1608 DO gamma = 1, 3
1609 blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1610 = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1611 matmul(achint(1:na, 1:np, i_1), transpose(bcint(1:nb, 1:np, delta + 1))) &
1612 *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1613 END DO
1614 END DO
1615 ELSE ! r_delta * V * (R^eta_gamma - R^nu_gamma)
1616 DO delta = 1, 3
1617 DO gamma = 1, 3
1618 blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1619 = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1620 matmul(achint(1:na, 1:np, delta + 1), transpose(bcint(1:nb, 1:np, i_1))) &
1621 *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1622 END DO
1623 END DO
1624 END IF
1625
1626!$ CALL omp_unset_lock(locks(hash))
1627 EXIT ! We have found a match and there can be only one single match
1628 END IF
1629 END DO
1630 END DO
1631 END DO
1632 DO delta = 1, 3
1633 DO gamma = 1, 3
1634 NULLIFY (blocks_rv(gamma, delta)%block)
1635 END DO
1636 END DO
1637 DEALLOCATE (blocks_rv)
1638 END DO
1639
1640!$OMP DO
1641!$ DO lock_num = 1, nlock
1642!$ call omp_destroy_lock(locks(lock_num))
1643!$ END DO
1644!$OMP END DO
1645
1646!$OMP SINGLE
1647!$ DEALLOCATE (locks)
1648!$OMP END SINGLE NOWAIT
1649
1650!$OMP END PARALLEL
1651
1652 CALL release_sap_int(sap_int)
1653
1654 DEALLOCATE (basis_set)
1655
1656 CALL timestop(handle)
1657
1658 END SUBROUTINE build_com_vnl_giao
1659
1660END MODULE commutator_rpnl
basis_set_types
Definition basis_set_types.F:15
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
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition commutator_rpnl.F:14
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition commutator_rpnl.F:76
commutator_rpnl::build_com_vnl_giao
subroutine, public build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_rv, ref_point, cell, direction_or)
Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs.
Definition commutator_rpnl.F:1438
commutator_rpnl::build_com_nl_mag
subroutine, public build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
Definition commutator_rpnl.F:1198
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
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
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
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_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
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
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
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_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
sap_kind_types::alist_type
Definition sap_kind_types.F:58
sap_kind_types::sap_int_type
Definition sap_kind_types.F:64