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 80 : CALL timeset(routineN, handle)
127 :
128 80 : my_r_doublecom = .FALSE.
129 80 : my_r_rxvr = .FALSE.
130 80 : my_rxvr_r = .FALSE.
131 80 : my_rxrv = .FALSE.
132 80 : my_rrv = .FALSE.
133 80 : my_rv = .FALSE.
134 80 : my_rvr = .FALSE.
135 80 : my_rrv_vrr = .FALSE.
136 80 : IF (PRESENT(matrix_r_doublecom)) my_r_doublecom = .TRUE.
137 80 : IF (PRESENT(matrix_r_rxvr)) my_r_rxvr = .TRUE.
138 80 : IF (PRESENT(matrix_rxvr_r)) my_rxvr_r = .TRUE.
139 80 : IF (PRESENT(matrix_rxrv)) my_rxrv = .TRUE.
140 80 : IF (PRESENT(matrix_rrv)) my_rrv = .TRUE.
141 80 : IF (PRESENT(matrix_rv)) my_rv = .TRUE.
142 80 : IF (PRESENT(matrix_rvr)) my_rvr = .TRUE.
143 80 : IF (PRESENT(matrix_rrv_vrr)) my_rrv_vrr = .TRUE.
144 80 : 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 80 : natom = SIZE(particle_set)
149 :
150 80 : IF (my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) THEN
151 46 : 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 80 : order = 1
156 80 : 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 20 : IF (my_r_doublecom) THEN
160 6 : CPASSERT(PRESENT(pseudoatom))
161 : END IF
162 :
163 200 : periodic = ANY(cell%perd > 0)
164 80 : my_ref = .FALSE.
165 80 : IF (PRESENT(ref_point)) THEN
166 52 : IF (.NOT. periodic) THEN
167 16 : rf = ref_point
168 16 : 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 36 : IF (order > 1 .OR. my_rxrv .OR. my_rvr) 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 80 : nkind = SIZE(qs_kind_set)
178 :
179 : !sap_int needs to be shared as multiple threads need to access this
180 80 : NULLIFY (sap_int)
181 752 : ALLOCATE (sap_int(nkind*nkind))
182 592 : DO i = 1, nkind*nkind
183 512 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
184 592 : sap_int(i)%nalist = 0
185 : END DO
186 :
187 80 : 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 16 : 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 80 : CALL sap_sort(sap_int)
197 :
198 432 : ALLOCATE (basis_set(nkind))
199 272 : DO ikind = 1, nkind
200 192 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
201 272 : IF (ASSOCIATED(orb_basis_set)) THEN
202 192 : 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 80 : 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 80 : !$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 80 : CALL release_sap_int(sap_int)
1144 :
1145 80 : DEALLOCATE (basis_set)
1146 :
1147 80 : CALL timestop(handle)
1148 :
1149 264 : 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 8 : 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 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1186 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1187 8 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1188 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1189 8 : DIMENSION(:) :: basis_set
1190 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1191 8 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1192 :
1193 : !$ INTEGER(kind=omp_lock_kind), &
1194 8 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1195 : !$ INTEGER :: lock_num, hash
1196 : !$ INTEGER, PARAMETER :: nlock = 501
1197 :
1198 8 : ppnl_present = ASSOCIATED(sap_ppnl)
1199 8 : IF (.NOT. ppnl_present) RETURN
1200 :
1201 8 : CALL timeset(routineN, handle)
1202 :
1203 8 : my_ref = .FALSE.
1204 8 : IF (PRESENT(refpoint)) THEN
1205 8 : my_ref = .TRUE.
1206 8 : CPASSERT(PRESENT(cell))
1207 : END IF
1208 :
1209 8 : natom = SIZE(particle_set)
1210 8 : nkind = SIZE(qs_kind_set)
1211 :
1212 : ! allocate integral storage
1213 8 : NULLIFY (sap_int)
1214 56 : ALLOCATE (sap_int(nkind*nkind))
1215 40 : DO ind = 1, nkind*nkind
1216 32 : NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1217 40 : sap_int(ind)%nalist = 0
1218 : END DO
1219 :
1220 : ! build integrals over GTO + projector functions, refpoint actually
1221 8 : 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 8 : IF (my_ref) THEN
1224 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE., refpoint=refpoint, &
1225 8 : particle_set=particle_set, cell=cell)
1226 : ELSE
1227 0 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE.)
1228 : END IF
1229 :
1230 8 : CALL sap_sort(sap_int)
1231 :
1232 : ! get access to basis sets
1233 40 : ALLOCATE (basis_set(nkind))
1234 24 : DO ikind = 1, nkind
1235 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1236 24 : IF (ASSOCIATED(orb_basis_set)) THEN
1237 16 : basis_set(ikind)%gto_basis_set => orb_basis_set
1238 : ELSE
1239 0 : 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 8 : !$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 8 : DEALLOCATE (basis_set)
1382 8 : CALL release_sap_int(sap_int)
1383 :
1384 8 : CALL timestop(handle)
1385 :
1386 16 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1428 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1429 0 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1430 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1431 0 : DIMENSION(:) :: basis_set
1432 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1433 0 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1434 :
1435 : !$ INTEGER(kind=omp_lock_kind), &
1436 0 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1437 : !$ INTEGER :: lock_num, hash
1438 : !$ INTEGER, PARAMETER :: nlock = 501
1439 :
1440 0 : ppnl_present = ASSOCIATED(sap_ppnl)
1441 0 : IF (.NOT. ppnl_present) RETURN
1442 :
1443 0 : CALL timeset(routineN, handle)
1444 :
1445 0 : natom = SIZE(particle_set)
1446 :
1447 0 : my_ref = .FALSE.
1448 0 : IF (PRESENT(ref_point)) THEN
1449 0 : CPASSERT(PRESENT(cell)) ! need cell as well if refpoint is provided
1450 0 : rf = ref_point
1451 0 : my_ref = .TRUE.
1452 : END IF
1453 :
1454 0 : nkind = SIZE(qs_kind_set)
1455 :
1456 : ! sap_int needs to be shared as multiple threads need to access this
1457 0 : NULLIFY (sap_int)
1458 0 : ALLOCATE (sap_int(nkind*nkind))
1459 0 : DO i = 1, nkind*nkind
1460 0 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1461 0 : sap_int(i)%nalist = 0
1462 : END DO
1463 :
1464 0 : order = 1
1465 0 : 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 0 : particle_set=particle_set, cell=cell)
1469 : ELSE
1470 0 : 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 0 : CALL sap_sort(sap_int)
1475 :
1476 0 : ALLOCATE (basis_set(nkind))
1477 0 : DO ikind = 1, nkind
1478 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1479 0 : IF (ASSOCIATED(orb_basis_set)) THEN
1480 0 : basis_set(ikind)%gto_basis_set => orb_basis_set
1481 : ELSE
1482 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1483 : END IF
1484 : END DO
1485 :
1486 0 : 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 0 : !$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 0 : CALL release_sap_int(sap_int)
1618 :
1619 0 : DEALLOCATE (basis_set)
1620 :
1621 0 : CALL timestop(handle)
1622 :
1623 0 : END SUBROUTINE build_com_vnl_giao
1624 :
1625 : END MODULE commutator_rpnl
|
