Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
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 : ! **************************************************************************************************
14 : MODULE 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 :
48 : CONTAINS
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 208 : SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, &
75 312 : 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 104 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
108 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
109 104 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_rrv, blocks_rrv_vrr, blocks_rv, &
110 104 : blocks_rvr, blocks_rxrv
111 104 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_r_doublecom, blocks_r_rxvr, &
112 104 : blocks_rxvr_r
113 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
114 104 : DIMENSION(:) :: basis_set
115 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
116 104 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
117 :
118 : !$ INTEGER(kind=omp_lock_kind), &
119 104 : !$ ALLOCATABLE, DIMENSION(:) :: locks
120 : !$ INTEGER :: lock_num, hash
121 : !$ INTEGER, PARAMETER :: nlock = 501
122 :
123 104 : ppnl_present = ASSOCIATED(sap_ppnl)
124 104 : IF (.NOT. ppnl_present) RETURN
125 :
126 82 : CALL timeset(routineN, handle)
127 :
128 82 : my_r_doublecom = .FALSE.
129 82 : my_r_rxvr = .FALSE.
130 82 : my_rxvr_r = .FALSE.
131 82 : my_rxrv = .FALSE.
132 82 : my_rrv = .FALSE.
133 82 : my_rv = .FALSE.
134 82 : my_rvr = .FALSE.
135 82 : my_rrv_vrr = .FALSE.
136 82 : IF (PRESENT(matrix_r_doublecom)) my_r_doublecom = .TRUE.
137 82 : IF (PRESENT(matrix_r_rxvr)) my_r_rxvr = .TRUE.
138 82 : IF (PRESENT(matrix_rxvr_r)) my_rxvr_r = .TRUE.
139 82 : IF (PRESENT(matrix_rxrv)) my_rxrv = .TRUE.
140 82 : IF (PRESENT(matrix_rrv)) my_rrv = .TRUE.
141 82 : IF (PRESENT(matrix_rv)) my_rv = .TRUE.
142 82 : IF (PRESENT(matrix_rvr)) my_rvr = .TRUE.
143 82 : IF (PRESENT(matrix_rrv_vrr)) my_rrv_vrr = .TRUE.
144 82 : 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 0 : CPABORT('No dbcsr matrix provided for commutator calculation!')
146 : END IF
147 :
148 82 : natom = SIZE(particle_set)
149 :
150 82 : IF (my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) THEN
151 48 : order = 2
152 48 : CPASSERT(PRESENT(ref_point)) ! need reference point for r x [r,Vnl] and [rr,Vnl]
153 34 : ELSE IF (my_rvr .OR. my_rrv_vrr) THEN
154 6 : order = 2
155 : ELSE
156 28 : order = 1
157 : END IF
158 :
159 : ! When we want the double commutator [[Vnl, r], r], we also want to fix the pseudoatom
160 82 : IF (my_r_doublecom) THEN
161 6 : CPASSERT(PRESENT(pseudoatom))
162 : END IF
163 :
164 208 : periodic = ANY(cell%perd > 0)
165 82 : my_ref = .FALSE.
166 82 : IF (PRESENT(ref_point)) THEN
167 54 : IF (.NOT. periodic) THEN
168 18 : rf = ref_point
169 18 : my_ref = .TRUE.
170 : ELSE ! use my_ref = False in periodic case, corresponds to distributed ref point
171 36 : IF (order > 1) THEN
172 36 : 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 82 : nkind = SIZE(qs_kind_set)
178 :
179 : !sap_int needs to be shared as multiple threads need to access this
180 82 : NULLIFY (sap_int)
181 766 : ALLOCATE (sap_int(nkind*nkind))
182 602 : DO i = 1, nkind*nkind
183 520 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
184 602 : sap_int(i)%nalist = 0
185 : END DO
186 :
187 82 : 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 18 : particle_set=particle_set, cell=cell)
191 : ELSE
192 64 : 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 82 : CALL sap_sort(sap_int)
197 :
198 442 : ALLOCATE (basis_set(nkind))
199 278 : DO ikind = 1, nkind
200 196 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
201 278 : IF (ASSOCIATED(orb_basis_set)) THEN
202 196 : basis_set(ikind)%gto_basis_set => orb_basis_set
203 : ELSE
204 0 : 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 82 : 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 82 : !$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 82 : CALL release_sap_int(sap_int)
1179 :
1180 82 : DEALLOCATE (basis_set)
1181 :
1182 82 : CALL timestop(handle)
1183 :
1184 268 : 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 8 : 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 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1221 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1222 8 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1223 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1224 8 : DIMENSION(:) :: basis_set
1225 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1226 8 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1227 :
1228 : !$ INTEGER(kind=omp_lock_kind), &
1229 8 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1230 : !$ INTEGER :: lock_num, hash
1231 : !$ INTEGER, PARAMETER :: nlock = 501
1232 :
1233 8 : ppnl_present = ASSOCIATED(sap_ppnl)
1234 8 : IF (.NOT. ppnl_present) RETURN
1235 :
1236 8 : CALL timeset(routineN, handle)
1237 :
1238 8 : my_ref = .FALSE.
1239 8 : IF (PRESENT(refpoint)) THEN
1240 8 : my_ref = .TRUE.
1241 8 : CPASSERT(PRESENT(cell))
1242 : END IF
1243 :
1244 8 : natom = SIZE(particle_set)
1245 8 : nkind = SIZE(qs_kind_set)
1246 :
1247 : ! allocate integral storage
1248 8 : NULLIFY (sap_int)
1249 56 : ALLOCATE (sap_int(nkind*nkind))
1250 40 : DO ind = 1, nkind*nkind
1251 32 : NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1252 40 : sap_int(ind)%nalist = 0
1253 : END DO
1254 :
1255 : ! build integrals over GTO + projector functions, refpoint actually
1256 8 : 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 8 : IF (my_ref) THEN
1259 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE., refpoint=refpoint, &
1260 8 : particle_set=particle_set, cell=cell)
1261 : ELSE
1262 0 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE.)
1263 : END IF
1264 :
1265 8 : CALL sap_sort(sap_int)
1266 :
1267 : ! get access to basis sets
1268 40 : ALLOCATE (basis_set(nkind))
1269 24 : DO ikind = 1, nkind
1270 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1271 24 : IF (ASSOCIATED(orb_basis_set)) THEN
1272 16 : basis_set(ikind)%gto_basis_set => orb_basis_set
1273 : ELSE
1274 0 : 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 8 : !$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 8 : DEALLOCATE (basis_set)
1417 8 : CALL release_sap_int(sap_int)
1418 :
1419 8 : CALL timestop(handle)
1420 :
1421 16 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1463 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1464 0 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1465 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1466 0 : DIMENSION(:) :: basis_set
1467 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1468 0 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1469 :
1470 : !$ INTEGER(kind=omp_lock_kind), &
1471 0 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1472 : !$ INTEGER :: lock_num, hash
1473 : !$ INTEGER, PARAMETER :: nlock = 501
1474 :
1475 0 : ppnl_present = ASSOCIATED(sap_ppnl)
1476 0 : IF (.NOT. ppnl_present) RETURN
1477 :
1478 0 : CALL timeset(routineN, handle)
1479 :
1480 0 : natom = SIZE(particle_set)
1481 :
1482 0 : my_ref = .FALSE.
1483 0 : IF (PRESENT(ref_point)) THEN
1484 0 : CPASSERT(PRESENT(cell)) ! need cell as well if refpoint is provided
1485 0 : rf = ref_point
1486 0 : my_ref = .TRUE.
1487 : END IF
1488 :
1489 0 : nkind = SIZE(qs_kind_set)
1490 :
1491 : ! sap_int needs to be shared as multiple threads need to access this
1492 0 : NULLIFY (sap_int)
1493 0 : ALLOCATE (sap_int(nkind*nkind))
1494 0 : DO i = 1, nkind*nkind
1495 0 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1496 0 : sap_int(i)%nalist = 0
1497 : END DO
1498 :
1499 0 : order = 1
1500 0 : 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 0 : particle_set=particle_set, cell=cell)
1504 : ELSE
1505 0 : 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 0 : CALL sap_sort(sap_int)
1510 :
1511 0 : ALLOCATE (basis_set(nkind))
1512 0 : DO ikind = 1, nkind
1513 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1514 0 : IF (ASSOCIATED(orb_basis_set)) THEN
1515 0 : basis_set(ikind)%gto_basis_set => orb_basis_set
1516 : ELSE
1517 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1518 : END IF
1519 : END DO
1520 :
1521 0 : 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 0 : !$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 0 : CALL release_sap_int(sap_int)
1653 :
1654 0 : DEALLOCATE (basis_set)
1655 :
1656 0 : CALL timestop(handle)
1657 :
1658 0 : END SUBROUTINE build_com_vnl_giao
1659 :
1660 : END MODULE commutator_rpnl
|
