(git:db0abe2)
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 cpassert(PRESENT(ref_point)) ! need reference point for r x [r,Vnl] and [rr,Vnl]
152 END IF
153
154 ! Only products with two coordinates on the same side need second-order SAP moments.
155 order = 1
156 IF (my_rrv .OR. my_rrv_vrr .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) order = 2
157
158 ! When we want the double commutator [[Vnl, r], r], we also want to fix the pseudoatom
159 IF (my_r_doublecom) THEN
160 cpassert(PRESENT(pseudoatom))
161 END IF
162
163 periodic = any(cell%perd > 0)
164 my_ref = .false.
165 IF (PRESENT(ref_point)) THEN
166 IF (.NOT. periodic) THEN
167 rf = ref_point
168 my_ref = .true.
169 ELSE ! use my_ref = False in periodic case, corresponds to distributed ref point
170 ! Reference-point diagnostics depend on the operator, not the SAP moment order.
171 IF (order > 1 .OR. my_rxrv .OR. my_rvr) 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 ! The same-side quadratic terms cancel in each cross-product component.
509 ! x-component (y [z,Vnl] - z [y, Vnl])
510 IF (iatom <= jatom) THEN
511 ! -yVz
512 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
513 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4)))
514 ! zVy
515 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
516 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 3)))
517 ELSE
518 ! -yVz
519 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
520 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4)))
521 ! zVy
522 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
523 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 3)))
524 END IF
525
526 ! y-component (z [x,Vnl] - x [z, Vnl])
527 IF (iatom <= jatom) THEN
528 ! -zVx
529 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
530 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 2)))
531 ! xVz
532 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
533 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4)))
534 ELSE
535 ! -zVx
536 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
537 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 2)))
538 ! xVz
539 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
540 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4)))
541 END IF
542
543 ! z-component (x [y,Vnl] - y [x, Vnl])
544 IF (iatom <= jatom) THEN
545 ! -xVy
546 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
547 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3)))
548 ! yVx
549 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
550 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 2)))
551 ELSE
552 ! -xVy
553 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
554 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3)))
555 ! yVx
556 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
557 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 2)))
558 END IF
559 END IF
560
561 IF (my_rrv) THEN
562 ! r_alpha * r_beta * Vnl
563 IF (iatom <= jatom) THEN
564 ! xxV
565 blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) + &
566 matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1)))
567 ! xyV
568 blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) + &
569 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
570 ! xzV
571 blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) + &
572 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
573 ! yyV
574 blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) + &
575 matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1)))
576 ! yzV
577 blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) + &
578 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
579 ! zzV
580 blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) + &
581 matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1)))
582 ELSE
583 ! xxV
584 blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) + &
585 matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1)))
586 ! xyV
587 blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) + &
588 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
589 ! xzV
590 blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) + &
591 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
592 ! yyV
593 blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) + &
594 matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1)))
595 ! yzV
596 blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) + &
597 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
598 ! zzV
599 blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) + &
600 matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1)))
601 END IF
602
603 ! - Vnl * r_alpha * r_beta
604 IF (iatom <= jatom) THEN
605 ! -Vxx
606 blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) - &
607 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5)))
608 ! -Vxy
609 blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) - &
610 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6)))
611 ! -Vxz
612 blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) - &
613 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7)))
614 ! -Vyy
615 blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) - &
616 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8)))
617 ! -Vyz
618 blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) - &
619 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9)))
620 ! -Vzz
621 blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) - &
622 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10)))
623 ELSE
624 ! -Vxx
625 blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) - &
626 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5)))
627 ! -Vxy
628 blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) - &
629 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6)))
630 ! -Vxz
631 blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) - &
632 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7)))
633 ! -Vyy
634 blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) - &
635 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8)))
636 ! -Vyz
637 blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) - &
638 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9)))
639 ! -Vzz
640 blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) - &
641 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10)))
642 END IF
643 END IF
644
645 IF (my_rvr) THEN
646 ! r_alpha * Vnl * r_beta
647 IF (iatom <= jatom) THEN
648 ! xVx
649 blocks_rvr(1)%block(1:na, 1:nb) = blocks_rvr(1)%block(1:na, 1:nb) + &
650 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 2)))
651 ! xVy
652 blocks_rvr(2)%block(1:na, 1:nb) = blocks_rvr(2)%block(1:na, 1:nb) + &
653 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3)))
654 ! xVz
655 blocks_rvr(3)%block(1:na, 1:nb) = blocks_rvr(3)%block(1:na, 1:nb) + &
656 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4)))
657 ! yVy
658 blocks_rvr(4)%block(1:na, 1:nb) = blocks_rvr(4)%block(1:na, 1:nb) + &
659 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 3)))
660 ! yVz
661 blocks_rvr(5)%block(1:na, 1:nb) = blocks_rvr(5)%block(1:na, 1:nb) + &
662 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4)))
663 ! zVz
664 blocks_rvr(6)%block(1:na, 1:nb) = blocks_rvr(6)%block(1:na, 1:nb) + &
665 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 4)))
666 ELSE
667 ! xVx
668 blocks_rvr(1)%block(1:nb, 1:na) = blocks_rvr(1)%block(1:nb, 1:na) + &
669 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 2)))
670 ! xVy
671 blocks_rvr(2)%block(1:nb, 1:na) = blocks_rvr(2)%block(1:nb, 1:na) + &
672 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3)))
673 ! xVz
674 blocks_rvr(3)%block(1:nb, 1:na) = blocks_rvr(3)%block(1:nb, 1:na) + &
675 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4)))
676 ! yVy
677 blocks_rvr(4)%block(1:nb, 1:na) = blocks_rvr(4)%block(1:nb, 1:na) + &
678 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 3)))
679 ! yVz
680 blocks_rvr(5)%block(1:nb, 1:na) = blocks_rvr(5)%block(1:nb, 1:na) + &
681 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4)))
682 ! zVz
683 blocks_rvr(6)%block(1:nb, 1:na) = blocks_rvr(6)%block(1:nb, 1:na) + &
684 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 4)))
685 END IF
686 END IF
687
688 IF (my_rrv_vrr) THEN
689 ! r_alpha * r_beta * Vnl
690 IF (iatom <= jatom) THEN
691 ! xxV
692 blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
693 matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1)))
694 ! xyV
695 blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
696 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
697 ! xzV
698 blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
699 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
700 ! yyV
701 blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
702 matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1)))
703 ! yzV
704 blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
705 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
706 ! zzV
707 blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
708 matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1)))
709 ELSE
710 ! xxV
711 blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
712 matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1)))
713 ! xyV
714 blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
715 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
716 ! xzV
717 blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
718 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
719 ! yyV
720 blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
721 matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1)))
722 ! yzV
723 blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
724 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
725 ! zzV
726 blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
727 matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1)))
728 END IF
729 ! + Vnl * r_alpha * r_beta
730 IF (iatom <= jatom) THEN
731 ! +Vxx
732 blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
733 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5)))
734 ! +Vxy
735 blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
736 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6)))
737 ! +Vxz
738 blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
739 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7)))
740 ! +Vyy
741 blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
742 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8)))
743 ! +Vyz
744 blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
745 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9)))
746 ! +Vzz
747 blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
748 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10)))
749 ELSE
750 ! +Vxx
751 blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
752 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5)))
753 ! +Vxy
754 blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
755 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6)))
756 ! +Vxz
757 blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
758 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7)))
759 ! +Vyy
760 blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
761 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8)))
762 ! +Vyz
763 blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
764 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9)))
765 ! +Vzz
766 blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
767 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10)))
768 END IF
769 END IF
770
771 ! The indices are stored in i_1, i_x, ..., i_zzz
772
773 ! TODO: is this set to zero before?
774 IF (my_r_rxvr) THEN
775 ! beta = 1
776 ! matrix_r_rxvr(x, x) = x * y * V_nl * z - x * z * V_nl * y
777 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
778 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) + &
779 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
780 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
781 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) - &
782 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
783
784 ! matrix_r_rxvr(y, x) = x * z * V_nl * x - x * x * V_nl * z
785 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
786 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) + &
787 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_x)))
788 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
789 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) - &
790 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
791
792 ! matrix_r_rxvr(z, x) = x * x * V_nl * y - x * y * V_nl * x
793 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
794 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) + &
795 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
796 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
797 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) - &
798 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_x)))
799
800 ! beta = 2
801 ! matrix_r_rxvr(x, y) = y * y * V_nl * z - y * z * V_nl * y
802 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
803 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) + &
804 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
805 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
806 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) - &
807 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_y)))
808
809 ! matrix_r_rxvr(y, y) = y * z * V_nl * x - y * x * V_nl * z
810 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
811 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) + &
812 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
813 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
814 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) - &
815 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
816
817 ! matrix_r_rxvr(z, y) = y * x * V_nl * y - y * y * V_nl * x
818 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
819 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) + &
820 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_y)))
821 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
822 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) - &
823 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
824
825 ! beta = 3
826 ! matrix_r_rxvr(x, z) = z * y * V_nl * z - z * z * V_nl * y
827 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
828 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) + &
829 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_z)))
830 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
831 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) - &
832 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
833
834 ! matrix_r_rxvr(y, z) = z * z * V_nl * x - z * x * V_nl * z
835 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
836 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) + &
837 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
838 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
839 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) - &
840 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_z)))
841
842 ! matrix_r_rxvr(z, z) = z * x * V_nl * y - z * y * V_nl * x
843 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
844 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) + &
845 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
846 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
847 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) - &
848 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
849
850 END IF ! my_r_rxvr
851
852 ! The indices are stored in i_1, i_x, ..., i_zzz
853 ! This will put into blocks_rxvr_r
854 ! matrix_rxvr_r(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
855 ! r_gamma * V_nl * r_delta * r_beta
856 IF (my_rxvr_r) THEN
857 ! beta = 1
858 ! matrix_rxvr_r(x, x) = yV zx - zV yx
859 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
860 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) + &
861 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
862 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
863 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) - &
864 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
865
866 ! matrix_rxvr_r(y, x) = zV xx - xV zx
867 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
868 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) + &
869 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
870 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
871 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) - &
872 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zx)))
873
874 ! matrix_rxvr_r(z, x) = xV yx - yV xx
875 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
876 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) + &
877 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yx)))
878 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
879 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) - &
880 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
881
882 ! beta = 2
883 ! matrix_rxvr_r(x, y) = yV zy - zV yy
884 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
885 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) + &
886 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zy)))
887 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
888 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) - &
889 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
890
891 ! matrix_rxvr_r(y, y) = zV xy - xV zy
892 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
893 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) + &
894 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
895 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
896 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) - &
897 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
898
899 ! matrix_rxvr_r(z, y) = xV yy - yV xy
900 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
901 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) + &
902 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
903 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
904 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) - &
905 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xy)))
906
907 ! beta = 3
908 ! matrix_rxvr_r(x, z) = yV zz - zV yz
909 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
910 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) + &
911 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
912 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
913 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) - &
914 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yz)))
915
916 ! matrix_rxvr_r(y, z) = zV xz - xV zz
917 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
918 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) + &
919 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xz)))
920 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
921 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) - &
922 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
923
924 ! matrix_rxvr_r(z, z) = xV yz - yV xz
925 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
926 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) + &
927 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
928 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
929 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) - &
930 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
931
932 END IF ! my_rxvr_r
933
934 ! matrix_r_doublecom(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
935 ! gamma V^pseudoatom beta delta - gamma beta V^pseudoatom delta
936
937 IF (my_r_doublecom) THEN
938 ! beta = 1
939 ! matrix_r_doublecom(x, x) = yV xz - zV xy - yxV z + zxV y
940 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
941 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
942 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
943 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
944 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
945 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
946 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
947 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
948 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
949 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
950 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
951 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
952
953 ! matrix_r_doublecom(y, x) = zV xx - xV xz - zxV x + xxV z
954 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
955 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
956 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
957 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
958 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
959 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xz)))
960 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
961 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
962 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_x)))
963 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
964 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
965 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
966
967 ! matrix_r_doublecom(z, x) = xV xy - yV xx - xxV y + yxV x
968 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
969 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
970 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xy)))
971 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
972 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
973 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
974 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
975 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
976 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
977 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
978 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
979 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_x)))
980
981 ! beta = 2
982 ! matrix_r_doublecom(x, y) = yV yz - zV yy - yyV z + zyV y
983 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
984 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
985 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yz)))
986 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
987 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
988 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
989 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
990 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
991 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
992 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
993 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
994 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_y)))
995
996 ! matrix_r_doublecom(y, y) = zV yx - xV yz - zyV x + xyV z
997 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
998 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
999 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
1000 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1001 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1002 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
1003 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1004 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1005 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
1006 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1007 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1008 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
1009
1010 ! matrix_r_doublecom(z, y) = xV yy - yV yx - xyV y + yyV x
1011 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1012 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1013 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
1014 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1015 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1016 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yx)))
1017 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1018 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1019 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_y)))
1020 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1021 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1022 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
1023
1024 ! beta = 3
1025 ! matrix_r_doublecom(x, z) = yV zz - zV zy - yzV z + zzV y
1026 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1027 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1028 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
1029 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1030 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1031 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zy)))
1032 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1033 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1034 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_z)))
1035 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1036 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1037 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
1038
1039 ! matrix_r_doublecom(y, z) = zV zx - xV zz - zzV x + xzV z
1040 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1041 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1042 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zx)))
1043 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1044 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1045 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
1046 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1047 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1048 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
1049 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1050 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1051 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_z)))
1052
1053 ! matrix_r_doublecom(z, z) = xV zy - yV zx - xzV y + yzV x
1054 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1055 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1056 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
1057 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1058 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1059 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
1060 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1061 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1062 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
1063 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1064 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1065 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
1066
1067 END IF ! my_r_doublecom
1068!$ CALL omp_unset_lock(locks(hash))
1069 EXIT ! We have found a match and there can be only one single match
1070 END IF
1071 END DO
1072 END DO
1073 END DO
1074 END IF
1075 IF (my_rv) THEN
1076 DO ind = 1, 3
1077 NULLIFY (blocks_rv(ind)%block)
1078 END DO
1079 DEALLOCATE (blocks_rv)
1080 END IF
1081 IF (my_rxrv) THEN
1082 DO ind = 1, 3
1083 NULLIFY (blocks_rxrv(ind)%block)
1084 END DO
1085 DEALLOCATE (blocks_rxrv)
1086 END IF
1087 IF (my_rrv) THEN
1088 DO ind = 1, 6
1089 NULLIFY (blocks_rrv(ind)%block)
1090 END DO
1091 DEALLOCATE (blocks_rrv)
1092 END IF
1093 IF (my_rvr) THEN
1094 DO ind = 1, 6
1095 NULLIFY (blocks_rvr(ind)%block)
1096 END DO
1097 DEALLOCATE (blocks_rvr)
1098 END IF
1099 IF (my_rrv_vrr) THEN
1100 DO ind = 1, 6
1101 NULLIFY (blocks_rrv_vrr(ind)%block)
1102 END DO
1103 DEALLOCATE (blocks_rrv_vrr)
1104 END IF
1105 IF (my_r_rxvr) THEN
1106 DO ind = 1, 3
1107 DO ind2 = 1, 3
1108 NULLIFY (blocks_r_rxvr(ind, ind2)%block)
1109 END DO
1110 END DO
1111 DEALLOCATE (blocks_r_rxvr)
1112 END IF
1113 IF (my_rxvr_r) THEN
1114 DO ind = 1, 3
1115 DO ind2 = 1, 3
1116 NULLIFY (blocks_rxvr_r(ind, ind2)%block)
1117 END DO
1118 END DO
1119 DEALLOCATE (blocks_rxvr_r)
1120 END IF
1121 IF (my_r_doublecom) THEN
1122 DO ind = 1, 3
1123 DO ind2 = 1, 3
1124 NULLIFY (blocks_r_doublecom(ind, ind2)%block)
1125 END DO
1126 END DO
1127 DEALLOCATE (blocks_r_doublecom)
1128 END IF
1129 END DO
1130
1131!$OMP DO
1132!$ DO lock_num = 1, nlock
1133!$ call omp_destroy_lock(locks(lock_num))
1134!$ END DO
1135!$OMP END DO
1136
1137!$OMP SINGLE
1138!$ DEALLOCATE (locks)
1139!$OMP END SINGLE NOWAIT
1140
1141!$OMP END PARALLEL
1142
1143 CALL release_sap_int(sap_int)
1144
1145 DEALLOCATE (basis_set)
1146
1147 CALL timestop(handle)
1148
1149 END SUBROUTINE build_com_mom_nl
1150
1151! **************************************************************************************************
1152!> \brief calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
1153!> \param qs_kind_set ...
1154!> \param sab_all ...
1155!> \param sap_ppnl ...
1156!> \param eps_ppnl ...
1157!> \param particle_set ...
1158!> \param matrix_mag_nl ...
1159!> \param refpoint ...
1160!> \param cell ...
1161! **************************************************************************************************
1162 SUBROUTINE build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
1163
1164 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1165 POINTER :: qs_kind_set
1166 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1167 INTENT(IN), POINTER :: sab_all, sap_ppnl
1168 REAL(kind=dp), INTENT(IN) :: eps_ppnl
1169 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1170 POINTER :: particle_set
1171 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1172 POINTER :: matrix_mag_nl
1173 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
1174 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1175
1176 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_nl_mag'
1177
1178 INTEGER :: handle, iab, iac, iatom, ibc, icol, &
1179 ikind, ind, irow, jatom, jkind, kac, &
1180 kbc, kkind, na, natom, nb, nkind, np, &
1181 order, slot
1182 INTEGER, DIMENSION(3) :: cell_b
1183 LOGICAL :: found, go, my_ref, ppnl_present
1184 REAL(kind=dp), DIMENSION(3) :: r_b, r_ps, rab
1185 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1186 TYPE(alist_type), POINTER :: alist_ac, alist_bc
1187 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1188 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1189 DIMENSION(:) :: basis_set
1190 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1191 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1192
1193!$ INTEGER(kind=omp_lock_kind), &
1194!$ ALLOCATABLE, DIMENSION(:) :: locks
1195!$ INTEGER :: lock_num, hash
1196!$ INTEGER, PARAMETER :: nlock = 501
1197
1198 ppnl_present = ASSOCIATED(sap_ppnl)
1199 IF (.NOT. ppnl_present) RETURN
1200
1201 CALL timeset(routinen, handle)
1202
1203 my_ref = .false.
1204 IF (PRESENT(refpoint)) THEN
1205 my_ref = .true.
1206 cpassert(PRESENT(cell))
1207 END IF
1208
1209 natom = SIZE(particle_set)
1210 nkind = SIZE(qs_kind_set)
1211
1212 ! allocate integral storage
1213 NULLIFY (sap_int)
1214 ALLOCATE (sap_int(nkind*nkind))
1215 DO ind = 1, nkind*nkind
1216 NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1217 sap_int(ind)%nalist = 0
1218 END DO
1219
1220 ! build integrals over GTO + projector functions, refpoint actually
1221 order = 1 ! only need first moments (x, y, z)
1222 ! refpoint actually does not matter in this case, i. e. (order = 1 .and. commutator)
1223 IF (my_ref) THEN
1224 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=refpoint, &
1225 particle_set=particle_set, cell=cell)
1226 ELSE
1227 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1228 END IF
1229
1230 CALL sap_sort(sap_int)
1231
1232 ! get access to basis sets
1233 ALLOCATE (basis_set(nkind))
1234 DO ikind = 1, nkind
1235 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1236 IF (ASSOCIATED(orb_basis_set)) THEN
1237 basis_set(ikind)%gto_basis_set => orb_basis_set
1238 ELSE
1239 NULLIFY (basis_set(ikind)%gto_basis_set)
1240 END IF
1241 END DO
1242
1243!$OMP PARALLEL &
1244!$OMP DEFAULT (NONE) &
1245!$OMP SHARED (basis_set, matrix_mag_nl, sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
1246!$OMP particle_set, my_ref, refpoint) &
1247!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, lock_num, &
1248!$OMP iab, irow, icol, blocks_mag, r_ps, r_b, go, hash, &
1249!$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
1250!$OMP achint, bchint, na, np, nb, kkind, kac, kbc)
1251
1252!$OMP SINGLE
1253!$ ALLOCATE (locks(nlock))
1254!$OMP END SINGLE
1255
1256!$OMP DO
1257!$ DO lock_num = 1, nlock
1258!$ call omp_init_lock(locks(lock_num))
1259!$ END DO
1260!$OMP END DO
1261
1262!$OMP DO SCHEDULE(GUIDED)
1263 DO slot = 1, sab_all(1)%nl_size
1264 ! get indices
1265 ikind = sab_all(1)%nlist_task(slot)%ikind
1266 jkind = sab_all(1)%nlist_task(slot)%jkind
1267 iatom = sab_all(1)%nlist_task(slot)%iatom
1268 jatom = sab_all(1)%nlist_task(slot)%jatom
1269 cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1270 rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1271
1272 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1273 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1274 iab = ikind + nkind*(jkind - 1)
1275
1276 IF (iatom <= jatom) THEN
1277 irow = iatom
1278 icol = jatom
1279 ELSE
1280 irow = jatom
1281 icol = iatom
1282 END IF
1283
1284 ! get blocks
1285 ALLOCATE (blocks_mag(3))
1286 DO ind = 1, 3
1287 CALL dbcsr_get_block_p(matrix_mag_nl(ind)%matrix, irow, icol, blocks_mag(ind)%block, found)
1288 END DO
1289
1290 go = (ASSOCIATED(blocks_mag(1)%block) .AND. ASSOCIATED(blocks_mag(2)%block) .AND. ASSOCIATED(blocks_mag(3)%block))
1291
1292 IF (go) THEN
1293 DO kkind = 1, nkind
1294 iac = ikind + nkind*(kkind - 1)
1295 ibc = jkind + nkind*(kkind - 1)
1296 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1297 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1298 CALL get_alist(sap_int(iac), alist_ac, iatom)
1299 CALL get_alist(sap_int(ibc), alist_bc, jatom)
1300 IF (.NOT. ASSOCIATED(alist_ac)) cycle
1301 IF (.NOT. ASSOCIATED(alist_bc)) cycle
1302 DO kac = 1, alist_ac%nclist
1303 DO kbc = 1, alist_bc%nclist
1304 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1305 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1306 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1307
1308 acint => alist_ac%clist(kac)%acint
1309 bcint => alist_bc%clist(kbc)%acint
1310 achint => alist_ac%clist(kac)%achint
1311 bchint => alist_bc%clist(kbc)%achint
1312 na = SIZE(acint, 1)
1313 np = SIZE(acint, 2)
1314 nb = SIZE(bcint, 1)
1315 ! Position of the pseudized atom
1316 r_ps = particle_set(alist_ac%clist(kac)%catom)%r
1317 r_b = refpoint
1318
1319!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1320!$ CALL omp_set_lock(locks(hash))
1321 ! assemble integrals
1322 IF (iatom <= jatom) THEN
1323 blocks_mag(1)%block(1:na, 1:nb) = blocks_mag(1)%block(1:na, 1:nb) + &
1324 (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1325 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) & ! R_y [V_nl, z]
1326 - (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1327 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) ! - R_z [V_nl, y]
1328 blocks_mag(2)%block(1:na, 1:nb) = blocks_mag(2)%block(1:na, 1:nb) + &
1329 (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1330 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) & ! R_z [V_nl, x]
1331 - (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1332 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) ! - R_x [V_nl, z]
1333 blocks_mag(3)%block(1:na, 1:nb) = blocks_mag(3)%block(1:na, 1:nb) + &
1334 (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1335 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) & ! R_x [V_nl, y]
1336 - (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1337 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) ! - R_y [V_nl, x]
1338 ELSE
1339 blocks_mag(1)%block(1:nb, 1:na) = blocks_mag(1)%block(1:nb, 1:na) + &
1340 (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1341 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) & ! R_y [V_nl, z]
1342 - (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1343 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) ! - R_z [V_nl, y]
1344 blocks_mag(2)%block(1:nb, 1:na) = blocks_mag(2)%block(1:nb, 1:na) + &
1345 (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1346 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) & ! R_z [V_nl, x]
1347 - (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1348 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) ! - R_x [V_nl, z]
1349 blocks_mag(3)%block(1:nb, 1:na) = blocks_mag(3)%block(1:nb, 1:na) + &
1350 (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1351 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) & ! R_x [V_nl, y]
1352 - (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1353 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) ! - R_y [V_nl, x]
1354 END IF
1355!$ CALL omp_unset_lock(locks(hash))
1356 EXIT ! We have found a match and there can be only one single match
1357 END IF
1358 END DO
1359 END DO
1360 END DO
1361 END IF
1362
1363 DO ind = 1, 3
1364 NULLIFY (blocks_mag(ind)%block)
1365 END DO
1366 DEALLOCATE (blocks_mag)
1367 END DO
1368
1369!$OMP DO
1370!$ DO lock_num = 1, nlock
1371!$ call omp_destroy_lock(locks(lock_num))
1372!$ END DO
1373!$OMP END DO
1374
1375!$OMP SINGLE
1376!$ DEALLOCATE (locks)
1377!$OMP END SINGLE NOWAIT
1378
1379!$OMP END PARALLEL
1380
1381 DEALLOCATE (basis_set)
1382 CALL release_sap_int(sap_int)
1383
1384 CALL timestop(handle)
1385
1386 END SUBROUTINE build_com_nl_mag
1387
1388! **************************************************************************************************
1389!> \brief Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs
1390!> \param qs_kind_set ...
1391!> \param sab_all ...
1392!> \param sap_ppnl ...
1393!> \param eps_ppnl ...
1394!> \param particle_set ...
1395!> \param matrix_rv ...
1396!> \param ref_point ...
1397!> \param cell ...
1398!> \param direction_Or If set to true: calculate Vnl * r_delta
1399!> Otherwise calculate r_delta * Vnl
1400! **************************************************************************************************
1401 SUBROUTINE build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, &
1402 matrix_rv, ref_point, cell, direction_Or)
1403
1404 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1405 POINTER :: qs_kind_set
1406 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1407 INTENT(IN), POINTER :: sab_all, sap_ppnl
1408 REAL(kind=dp), INTENT(IN) :: eps_ppnl
1409 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1410 POINTER :: particle_set
1411 TYPE(dbcsr_p_type), DIMENSION(:, :), &
1412 INTENT(INOUT), OPTIONAL, POINTER :: matrix_rv
1413 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
1414 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1415 LOGICAL :: direction_or
1416
1417 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_vnl_giao'
1418 INTEGER, PARAMETER :: i_1 = 1
1419
1420 INTEGER :: delta, gamma, handle, i, iab, iac, &
1421 iatom, ibc, icol, ikind, irow, j, &
1422 jatom, jkind, kac, kbc, kkind, na, &
1423 natom, nb, nkind, np, order, slot
1424 INTEGER, DIMENSION(3) :: cell_b
1425 LOGICAL :: found, my_ref, ppnl_present
1426 REAL(kind=dp), DIMENSION(3) :: rab, rf
1427 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1428 TYPE(alist_type), POINTER :: alist_ac, alist_bc
1429 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1430 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1431 DIMENSION(:) :: basis_set
1432 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1433 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1434
1435!$ INTEGER(kind=omp_lock_kind), &
1436!$ ALLOCATABLE, DIMENSION(:) :: locks
1437!$ INTEGER :: lock_num, hash
1438!$ INTEGER, PARAMETER :: nlock = 501
1439
1440 ppnl_present = ASSOCIATED(sap_ppnl)
1441 IF (.NOT. ppnl_present) RETURN
1442
1443 CALL timeset(routinen, handle)
1444
1445 natom = SIZE(particle_set)
1446
1447 my_ref = .false.
1448 IF (PRESENT(ref_point)) THEN
1449 cpassert(PRESENT(cell)) ! need cell as well if refpoint is provided
1450 rf = ref_point
1451 my_ref = .true.
1452 END IF
1453
1454 nkind = SIZE(qs_kind_set)
1455
1456 ! sap_int needs to be shared as multiple threads need to access this
1457 NULLIFY (sap_int)
1458 ALLOCATE (sap_int(nkind*nkind))
1459 DO i = 1, nkind*nkind
1460 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1461 sap_int(i)%nalist = 0
1462 END DO
1463
1464 order = 1
1465 IF (my_ref) THEN
1466 ! calculate integrals <a|x^n|p>
1467 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=rf, &
1468 particle_set=particle_set, cell=cell)
1469 ELSE
1470 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1471 END IF
1472
1473 ! *** Set up a sorting index
1474 CALL sap_sort(sap_int)
1475
1476 ALLOCATE (basis_set(nkind))
1477 DO ikind = 1, nkind
1478 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1479 IF (ASSOCIATED(orb_basis_set)) THEN
1480 basis_set(ikind)%gto_basis_set => orb_basis_set
1481 ELSE
1482 NULLIFY (basis_set(ikind)%gto_basis_set)
1483 END IF
1484 END DO
1485
1486 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_all)
1487 ! *** All integrals needed have been calculated and stored in sap_int
1488 ! *** We now calculate the commutator matrix elements
1489
1490!$OMP PARALLEL &
1491!$OMP DEFAULT (NONE) &
1492!$OMP SHARED (basis_set, matrix_rv, &
1493!$OMP sap_int, nkind, eps_ppnl, locks, sab_all, &
1494!$OMP particle_set, direction_Or) &
1495!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1496!$OMP iab, irow, icol, blocks_rv, &
1497!$OMP found, iac, ibc, alist_ac, alist_bc, &
1498!$OMP na, np, nb, kkind, kac, kbc, i, lock_num, &
1499!$OMP hash, natom, delta, gamma, achint, bchint, acint, bcint)
1500
1501!$OMP SINGLE
1502!$ ALLOCATE (locks(nlock))
1503!$OMP END SINGLE
1504
1505!$OMP DO
1506!$ DO lock_num = 1, nlock
1507!$ call omp_init_lock(locks(lock_num))
1508!$ END DO
1509!$OMP END DO
1510
1511!$OMP DO SCHEDULE(GUIDED)
1512
1513 DO slot = 1, sab_all(1)%nl_size
1514
1515 ikind = sab_all(1)%nlist_task(slot)%ikind
1516 jkind = sab_all(1)%nlist_task(slot)%jkind
1517 iatom = sab_all(1)%nlist_task(slot)%iatom
1518 jatom = sab_all(1)%nlist_task(slot)%jatom
1519 cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1520 rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1521
1522 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1523 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1524 iab = ikind + nkind*(jkind - 1)
1525
1526 irow = iatom
1527 icol = jatom
1528
1529 ! allocate blocks
1530 ALLOCATE (blocks_rv(3, 3))
1531
1532 ! get blocks
1533 DO i = 1, 3
1534 DO j = 1, 3
1535 CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, &
1536 blocks_rv(i, j)%block, found)
1537 blocks_rv(i, j)%block(:, :) = 0.0_dp
1538 cpassert(found)
1539 END DO
1540 END DO
1541
1542 ! loop over all kinds for projector atom
1543 ! < iatom | katom > h < katom | jatom >
1544 DO kkind = 1, nkind
1545 iac = ikind + nkind*(kkind - 1)
1546 ibc = jkind + nkind*(kkind - 1)
1547 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1548 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1549 CALL get_alist(sap_int(iac), alist_ac, iatom)
1550 CALL get_alist(sap_int(ibc), alist_bc, jatom)
1551 IF (.NOT. ASSOCIATED(alist_ac)) cycle
1552 IF (.NOT. ASSOCIATED(alist_bc)) cycle
1553 DO kac = 1, alist_ac%nclist
1554 DO kbc = 1, alist_bc%nclist
1555 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1556
1557 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1558 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1559 acint => alist_ac%clist(kac)%acint
1560 bcint => alist_bc%clist(kbc)%acint
1561 achint => alist_ac%clist(kac)%achint
1562 bchint => alist_bc%clist(kbc)%achint
1563 na = SIZE(acint, 1)
1564 np = SIZE(acint, 2)
1565 nb = SIZE(bcint, 1)
1566!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1567!$ CALL omp_set_lock(locks(hash))
1568
1569 !! The atom index is alist_ac%clist(kac)%catom
1570 ! The coordinate is particle_set(alist_ac%clist(kac)%catom)%r(:)
1571 IF (direction_or) THEN ! V * r_delta * (R^eta_gamma - R^nu_gamma)
1572 DO delta = 1, 3
1573 DO gamma = 1, 3
1574 blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1575 = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1576 matmul(achint(1:na, 1:np, i_1), transpose(bcint(1:nb, 1:np, delta + 1))) &
1577 *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1578 END DO
1579 END DO
1580 ELSE ! r_delta * V * (R^eta_gamma - R^nu_gamma)
1581 DO delta = 1, 3
1582 DO gamma = 1, 3
1583 blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1584 = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1585 matmul(achint(1:na, 1:np, delta + 1), transpose(bcint(1:nb, 1:np, i_1))) &
1586 *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1587 END DO
1588 END DO
1589 END IF
1590
1591!$ CALL omp_unset_lock(locks(hash))
1592 EXIT ! We have found a match and there can be only one single match
1593 END IF
1594 END DO
1595 END DO
1596 END DO
1597 DO delta = 1, 3
1598 DO gamma = 1, 3
1599 NULLIFY (blocks_rv(gamma, delta)%block)
1600 END DO
1601 END DO
1602 DEALLOCATE (blocks_rv)
1603 END DO
1604
1605!$OMP DO
1606!$ DO lock_num = 1, nlock
1607!$ call omp_destroy_lock(locks(lock_num))
1608!$ END DO
1609!$OMP END DO
1610
1611!$OMP SINGLE
1612!$ DEALLOCATE (locks)
1613!$OMP END SINGLE NOWAIT
1614
1615!$OMP END PARALLEL
1616
1617 CALL release_sap_int(sap_int)
1618
1619 DEALLOCATE (basis_set)
1620
1621 CALL timestop(handle)
1622
1623 END SUBROUTINE build_com_vnl_giao
1624
1625END MODULE commutator_rpnl
sap_kind_types::get_alist
Definition sap_kind_types.F:113
sap_kind_types::release_sap_int
Definition sap_kind_types.F:110
sap_kind_types::sap_sort
Definition sap_kind_types.F:116
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:1403
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:1163
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:483
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:425
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:61
sap_kind_types::sap_int_type
Definition sap_kind_types.F:67