(git:e7e05ae)
commutator_rpnl.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ai_moments, ONLY: moment
16  USE ai_overlap, ONLY: overlap
17  USE basis_set_types, ONLY: gto_basis_set_p_type,&
18  gto_basis_set_type
19  USE block_p_types, ONLY: block_p_type
20  USE cell_types, ONLY: cell_type
21  USE dbcsr_api, ONLY: dbcsr_get_block_p,&
22  dbcsr_p_type
23  USE external_potential_types, ONLY: gth_potential_p_type,&
24  gth_potential_type,&
25  sgp_potential_p_type,&
26  sgp_potential_type
27  USE kinds, ONLY: dp
28  USE orbital_pointers, ONLY: init_orbital_pointers,&
29  nco,&
30  ncoset
31  USE particle_types, ONLY: particle_type
32  USE qs_kind_types, ONLY: get_qs_kind,&
33  get_qs_kind_set,&
34  qs_kind_type
35  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
36  get_neighbor_list_set_p,&
37  neighbor_list_iterate,&
38  neighbor_list_iterator_create,&
39  neighbor_list_iterator_p_type,&
40  neighbor_list_iterator_release,&
41  neighbor_list_set_p_type
42  USE sap_kind_types, ONLY: alist_type,&
43  build_sap_ints,&
44  clist_type,&
45  get_alist,&
46  release_sap_int,&
47  sap_int_type,&
48  sap_sort
49 
50 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
51 !$ USE OMP_LIB, ONLY: omp_lock_kind, &
52 !$ omp_init_lock, omp_set_lock, &
53 !$ omp_unset_lock, omp_destroy_lock
54 
55 #include "./base/base_uses.f90"
56 
57  IMPLICIT NONE
58 
59  PRIVATE
60 
61  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'commutator_rpnl'
62 
63  PUBLIC :: build_com_rpnl, build_com_mom_nl, build_com_nl_mag, build_com_vnl_giao
64 
65 CONTAINS
66 
67 ! **************************************************************************************************
68 !> \brief ...
69 !> \param matrix_rv ...
70 !> \param qs_kind_set ...
71 !> \param sab_orb ...
72 !> \param sap_ppnl ...
73 !> \param eps_ppnl ...
74 ! **************************************************************************************************
75  SUBROUTINE build_com_rpnl(matrix_rv, qs_kind_set, sab_orb, sap_ppnl, eps_ppnl)
76 
77  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rv
78  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
79  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
80  POINTER :: sab_orb, sap_ppnl
81  REAL(kind=dp), INTENT(IN) :: eps_ppnl
82 
83  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_rpnl'
84 
85  INTEGER :: handle, i, iab, iac, iatom, ibc, icol, ikind, ilist, inode, irow, iset, jatom, &
86  jkind, jneighbor, kac, katom, kbc, kkind, l, lc_max, lc_min, ldai, ldsab, lppnl, maxco, &
87  maxder, maxl, maxlgto, maxlppnl, maxppnl, maxsgf, mepos, na, nb, ncoa, ncoc, nkind, &
88  nlist, nneighbor, nnode, np, nppnl, nprjc, nseta, nsgfa, nthread, prjc, sgfa
89  INTEGER, DIMENSION(3) :: cell_b, cell_c
90  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
91  nsgf_seta
92  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
93  LOGICAL :: found, gpot, ppnl_present, spot
94  REAL(kind=dp) :: dac, ppnl_radius
95  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work, sab, work
96  REAL(kind=dp), DIMENSION(1) :: rprjc, zetc
97  REAL(kind=dp), DIMENSION(3) :: rab, rac
98  REAL(kind=dp), DIMENSION(:), POINTER :: alpha_ppnl, set_radius_a
99  REAL(kind=dp), DIMENSION(:, :), POINTER :: cprj, rpgfa, sphi_a, vprj_ppnl, x_block, &
100  y_block, z_block, zeta
101  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
102  TYPE(alist_type), POINTER :: alist_ac, alist_bc
103  TYPE(clist_type), POINTER :: clist
104  TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
105  TYPE(gth_potential_type), POINTER :: gth_potential
106  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
107  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
108  TYPE(neighbor_list_iterator_p_type), &
109  DIMENSION(:), POINTER :: nl_iterator
110  TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
111  TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
112  TYPE(sgp_potential_type), POINTER :: sgp_potential
113 
114  CALL timeset(routinen, handle)
115 
116  ppnl_present = ASSOCIATED(sap_ppnl)
117 
118  IF (ppnl_present) THEN
119 
120  nkind = SIZE(qs_kind_set)
121 
122  CALL get_qs_kind_set(qs_kind_set, &
123  maxco=maxco, &
124  maxlgto=maxlgto, &
125  maxsgf=maxsgf, &
126  maxlppnl=maxlppnl, &
127  maxppnl=maxppnl)
128 
129  maxl = max(maxlgto, maxlppnl)
130  CALL init_orbital_pointers(maxl + 1)
131 
132  ldsab = max(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
133  ldai = ncoset(maxl + 1)
134 
135  !sap_int needs to be shared as multiple threads need to access this
136  ALLOCATE (sap_int(nkind*nkind))
137  DO i = 1, nkind*nkind
138  NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
139  sap_int(i)%nalist = 0
140  END DO
141 
142  !set up direct access to basis and potential
143  ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
144  DO ikind = 1, nkind
145  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
146  IF (ASSOCIATED(orb_basis_set)) THEN
147  basis_set(ikind)%gto_basis_set => orb_basis_set
148  ELSE
149  NULLIFY (basis_set(ikind)%gto_basis_set)
150  END IF
151  CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, &
152  sgp_potential=sgp_potential)
153  IF (ASSOCIATED(gth_potential)) THEN
154  gpotential(ikind)%gth_potential => gth_potential
155  NULLIFY (spotential(ikind)%sgp_potential)
156  ELSE IF (ASSOCIATED(sgp_potential)) THEN
157  spotential(ikind)%sgp_potential => sgp_potential
158  NULLIFY (gpotential(ikind)%gth_potential)
159  ELSE
160  NULLIFY (gpotential(ikind)%gth_potential)
161  NULLIFY (spotential(ikind)%sgp_potential)
162  END IF
163  END DO
164 
165  maxder = 4
166  nthread = 1
167 !$ nthread = omp_get_max_threads()
168 
169  !calculate the overlap integrals <a|p>
170  CALL neighbor_list_iterator_create(nl_iterator, sap_ppnl, nthread=nthread)
171 !$OMP PARALLEL &
172 !$OMP DEFAULT (NONE) &
173 !$OMP SHARED (nl_iterator, basis_set, spotential, gpotential, maxder, ncoset, &
174 !$OMP sap_int, nkind, ldsab, ldai, nco ) &
175 !$OMP PRIVATE (mepos, ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
176 !$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
177 !$OMP sphi_a, zeta, cprj, lppnl, nppnl, nprj_ppnl, &
178 !$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, &
179 !$OMP ncoc, rpgfa, vprj_ppnl, i, l, gpot, spot, &
180 !$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl)
181  mepos = 0
182 !$ mepos = omp_get_thread_num()
183 
184  ALLOCATE (sab(ldsab, ldsab, maxder), work(ldsab, ldsab, maxder))
185  sab = 0.0_dp
186  ALLOCATE (ai_work(ldai, ldai, 1))
187  ai_work = 0.0_dp
188 
189  DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
190  CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=kkind, iatom=iatom, &
191  jatom=katom, nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, cell=cell_c, r=rac)
192  iac = ikind + nkind*(kkind - 1)
193  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
194  gpot = ASSOCIATED(gpotential(kkind)%gth_potential)
195  spot = ASSOCIATED(spotential(kkind)%sgp_potential)
196  IF ((.NOT. gpot) .AND. (.NOT. spot)) cycle
197  ! get definition of basis set
198  first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf
199  la_max => basis_set(ikind)%gto_basis_set%lmax
200  la_min => basis_set(ikind)%gto_basis_set%lmin
201  npgfa => basis_set(ikind)%gto_basis_set%npgf
202  nseta = basis_set(ikind)%gto_basis_set%nset
203  nsgfa = basis_set(ikind)%gto_basis_set%nsgf
204  nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set
205  rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius
206  set_radius_a => basis_set(ikind)%gto_basis_set%set_radius
207  sphi_a => basis_set(ikind)%gto_basis_set%sphi
208  zeta => basis_set(ikind)%gto_basis_set%zet
209  nsgfa = basis_set(ikind)%gto_basis_set%nsgf
210 
211  ! get definition of PP projectors
212  IF (gpot) THEN
213  alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl
214  cprj => gpotential(kkind)%gth_potential%cprj
215  lppnl = gpotential(kkind)%gth_potential%lppnl
216  nppnl = gpotential(kkind)%gth_potential%nppnl
217  nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
218  ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
219  vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
220  ELSEIF (spot) THEN
221  cpabort('SGP not implemented')
222  ELSE
223  cpabort('PPNL unknown')
224  END IF
225 !$OMP CRITICAL(sap_int_critical)
226  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
227  sap_int(iac)%a_kind = ikind
228  sap_int(iac)%p_kind = kkind
229  sap_int(iac)%nalist = nlist
230  ALLOCATE (sap_int(iac)%alist(nlist))
231  DO i = 1, nlist
232  NULLIFY (sap_int(iac)%alist(i)%clist)
233  sap_int(iac)%alist(i)%aatom = 0
234  sap_int(iac)%alist(i)%nclist = 0
235  END DO
236  END IF
237  IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
238  sap_int(iac)%alist(ilist)%aatom = iatom
239  sap_int(iac)%alist(ilist)%nclist = nneighbor
240  ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
241  DO i = 1, nneighbor
242  sap_int(iac)%alist(ilist)%clist(i)%catom = 0
243  END DO
244  END IF
245 !$OMP END CRITICAL(sap_int_critical)
246  dac = sqrt(sum(rac*rac))
247  clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
248  clist%catom = katom
249  clist%cell = cell_c
250  clist%rac = rac
251  ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
252  clist%achint(nsgfa, nppnl, maxder))
253  clist%acint = 0._dp
254  clist%achint = 0._dp
255  clist%nsgf_cnt = 0
256  NULLIFY (clist%sgf_list)
257  DO iset = 1, nseta
258  ncoa = npgfa(iset)*ncoset(la_max(iset))
259  sgfa = first_sgfa(1, iset)
260  work = 0._dp
261  prjc = 1
262  DO l = 0, lppnl
263  nprjc = nprj_ppnl(l)*nco(l)
264  IF (nprjc == 0) cycle
265  rprjc(1) = ppnl_radius
266  IF (set_radius_a(iset) + rprjc(1) < dac) cycle
267  lc_max = l + 2*(nprj_ppnl(l) - 1)
268  lc_min = l
269  zetc(1) = alpha_ppnl(l)
270  ncoc = ncoset(lc_max)
271  ! Calculate the primitive overlap and dipole moment integrals
272  CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
273  lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab(:, :, 1), 0, .false., ai_work, ldai)
274  CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
275  lc_max, 1, zetc, rprjc, 1, rac, (/0._dp, 0._dp, 0._dp/), sab(:, :, 2:4))
276  ! *** Transformation step projector functions (cartesian->spherical) ***
277  DO i = 1, maxder
278  CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, 1, i), ldsab, &
279  cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, 1, i), ldsab)
280  END DO
281  prjc = prjc + nprjc
282  END DO
283  DO i = 1, maxder
284  ! Contraction step (basis functions)
285  CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
286  work(1, 1, i), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
287  ! Multiply with interaction matrix(h)
288  CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, &
289  vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa)
290  END DO
291  END DO
292  clist%maxac = maxval(abs(clist%acint(:, :, 1)))
293  clist%maxach = maxval(abs(clist%achint(:, :, 1)))
294  END DO
295 
296  DEALLOCATE (sab, ai_work, work)
297 !$OMP END PARALLEL
298  CALL neighbor_list_iterator_release(nl_iterator)
299 
300  ! *** Set up a sorting index
301  CALL sap_sort(sap_int)
302  ! *** All integrals needed have been calculated and stored in sap_int
303  ! *** We now calculate the Hamiltonian matrix elements
304  CALL neighbor_list_iterator_create(nl_iterator, sab_orb, nthread=nthread)
305 
306 !$OMP PARALLEL &
307 !$OMP DEFAULT (NONE) &
308 !$OMP SHARED (nl_iterator, basis_set, matrix_rv, &
309 !$OMP sap_int, nkind, eps_ppnl ) &
310 !$OMP PRIVATE (mepos, ikind, jkind, iatom, jatom, nlist, ilist, nnode, inode, cell_b, rab, &
311 !$OMP iab, irow, icol, x_block, y_block, z_block, &
312 !$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
313 !$OMP achint, bchint, na, np, nb, katom, rac, kkind, kac, kbc, i)
314 
315  mepos = 0
316 !$ mepos = omp_get_thread_num()
317 
318  DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
319  CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, &
320  jatom=jatom, nlist=nlist, ilist=ilist, nnode=nnode, inode=inode, cell=cell_b, r=rab)
321  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
322  IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
323  iab = ikind + nkind*(jkind - 1)
324 
325  ! *** Create matrix blocks for a new matrix block column ***
326  IF (iatom <= jatom) THEN
327  irow = iatom
328  icol = jatom
329  ELSE
330  irow = jatom
331  icol = iatom
332  END IF
333  CALL dbcsr_get_block_p(matrix_rv(1)%matrix, irow, icol, x_block, found)
334  CALL dbcsr_get_block_p(matrix_rv(2)%matrix, irow, icol, y_block, found)
335  CALL dbcsr_get_block_p(matrix_rv(3)%matrix, irow, icol, z_block, found)
336 
337  ! loop over all kinds for projector atom
338  IF (ASSOCIATED(x_block) .AND. ASSOCIATED(y_block) .AND. ASSOCIATED(z_block)) THEN
339  DO kkind = 1, nkind
340  iac = ikind + nkind*(kkind - 1)
341  ibc = jkind + nkind*(kkind - 1)
342  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
343  IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
344  CALL get_alist(sap_int(iac), alist_ac, iatom)
345  CALL get_alist(sap_int(ibc), alist_bc, jatom)
346  IF (.NOT. ASSOCIATED(alist_ac)) cycle
347  IF (.NOT. ASSOCIATED(alist_bc)) cycle
348  DO kac = 1, alist_ac%nclist
349  DO kbc = 1, alist_bc%nclist
350  IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
351  IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
352  IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
353  acint => alist_ac%clist(kac)%acint
354  bcint => alist_bc%clist(kbc)%acint
355  achint => alist_ac%clist(kac)%achint
356  bchint => alist_bc%clist(kbc)%achint
357  na = SIZE(acint, 1)
358  np = SIZE(acint, 2)
359  nb = SIZE(bcint, 1)
360 !$OMP CRITICAL(h_block_critical)
361  IF (iatom <= jatom) THEN
362  ! Vnl*r
363  CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 1), na, &
364  bcint(1, 1, 2), nb, 1.0_dp, x_block, SIZE(x_block, 1))
365  CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 1), na, &
366  bcint(1, 1, 3), nb, 1.0_dp, y_block, SIZE(y_block, 1))
367  CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 1), na, &
368  bcint(1, 1, 4), nb, 1.0_dp, z_block, SIZE(z_block, 1))
369  ! -r*Vnl
370  CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 2), na, &
371  bcint(1, 1, 1), nb, 1.0_dp, x_block, SIZE(x_block, 1))
372  CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 3), na, &
373  bcint(1, 1, 1), nb, 1.0_dp, y_block, SIZE(y_block, 1))
374  CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 4), na, &
375  bcint(1, 1, 1), nb, 1.0_dp, z_block, SIZE(z_block, 1))
376  ELSE
377  ! Vnl*r
378  CALL dgemm("N", "T", nb, na, np, 1.0_dp, bchint(1, 1, 2), nb, &
379  acint(1, 1, 1), na, 1.0_dp, x_block, SIZE(x_block, 1))
380  CALL dgemm("N", "T", nb, na, np, 1.0_dp, bchint(1, 1, 3), nb, &
381  acint(1, 1, 1), na, 1.0_dp, y_block, SIZE(y_block, 1))
382  CALL dgemm("N", "T", nb, na, np, 1.0_dp, bchint(1, 1, 4), nb, &
383  acint(1, 1, 1), na, 1.0_dp, z_block, SIZE(z_block, 1))
384  ! -r*Vnl
385  CALL dgemm("N", "T", nb, na, np, -1.0_dp, bchint(1, 1, 1), nb, &
386  acint(1, 1, 2), na, 1.0_dp, x_block, SIZE(x_block, 1))
387  CALL dgemm("N", "T", nb, na, np, -1.0_dp, bchint(1, 1, 1), nb, &
388  acint(1, 1, 3), na, 1.0_dp, y_block, SIZE(y_block, 1))
389  CALL dgemm("N", "T", nb, na, np, -1.0_dp, bchint(1, 1, 1), nb, &
390  acint(1, 1, 4), na, 1.0_dp, z_block, SIZE(z_block, 1))
391  END IF
392 !$OMP END CRITICAL(h_block_critical)
393  EXIT ! We have found a match and there can be only one single match
394  END IF
395  END DO
396  END DO
397  END DO
398  END IF
399  END DO
400 !$OMP END PARALLEL
401  CALL neighbor_list_iterator_release(nl_iterator)
402 
403  CALL release_sap_int(sap_int)
404 
405  DEALLOCATE (basis_set, gpotential, spotential)
406 
407  END IF !ppnl_present
408 
409  CALL timestop(handle)
410 
411  END SUBROUTINE build_com_rpnl
412 
413 ! **************************************************************************************************
414 !> \brief Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv)
415 !> or [rr,Vnl] (matrix_rrv) in AO basis.
416 !> Reference point is required for the two latter options
417 !> Update: Calculate rxVnlxr (matrix_rvr) and rxrxVnl + Vnlxrxr (matrix_rrv_vrr)
418 !> in AO basis. Added in the first place for current correction in
419 !> the VG formalism (first order wrt vector potential).
420 !> \param qs_kind_set ...
421 !> \param sab_all ...
422 !> \param sap_ppnl ...
423 !> \param eps_ppnl ...
424 !> \param particle_set ...
425 !> \param cell ...
426 !> \param matrix_rv ...
427 !> \param matrix_rxrv ...
428 !> \param matrix_rrv ...
429 !> \param matrix_rvr ...
430 !> \param matrix_rrv_vrr ...
431 !> \param matrix_r_rxvr ...
432 !> \param matrix_rxvr_r ...
433 !> \param matrix_r_doublecom ...
434 !> \param pseudoatom ...
435 !> \param ref_point ...
436 ! **************************************************************************************************
437  SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, &
438  matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
439 
440  TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
441  POINTER :: qs_kind_set
442  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
443  INTENT(IN), POINTER :: sab_all, sap_ppnl
444  REAL(kind=dp), INTENT(IN) :: eps_ppnl
445  TYPE(particle_type), DIMENSION(:), INTENT(IN), &
446  POINTER :: particle_set
447  TYPE(cell_type), INTENT(IN), POINTER :: cell
448  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
449  OPTIONAL :: matrix_rv, matrix_rxrv, matrix_rrv, &
450  matrix_rvr, matrix_rrv_vrr
451  TYPE(dbcsr_p_type), DIMENSION(:, :), &
452  INTENT(INOUT), OPTIONAL :: matrix_r_rxvr, matrix_rxvr_r, &
453  matrix_r_doublecom
454  INTEGER, INTENT(in), OPTIONAL :: pseudoatom
455  REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
456 
457  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_mom_nl'
458  INTEGER, PARAMETER :: i_x = 2, i_xx = 5, i_xy = 6, i_xz = 7, i_y = 3, i_yx = i_xy, i_yy = 8, &
459  i_yz = 9, i_z = 4, i_zx = i_xz, i_zy = i_yz, i_zz = 10
460 
461  INTEGER :: handle, i, iab, iac, iatom, ibc, icol, &
462  ikind, ind, ind2, irow, jatom, jkind, &
463  kac, kbc, kkind, na, natom, nb, nkind, &
464  np, order, slot
465  INTEGER, DIMENSION(3) :: cell_b
466  LOGICAL :: asso_r_doublecom, asso_r_rxvr, asso_rrv, asso_rrv_vrr, asso_rv, asso_rvr, &
467  asso_rxrv, asso_rxvr_r, do_symmetric, found, go, my_r_doublecom, my_r_rxvr, my_ref, &
468  my_rrv, my_rrv_vrr, my_rv, my_rvr, my_rxrv, my_rxvr_r, periodic, ppnl_present
469  REAL(kind=dp), DIMENSION(3) :: rab, rf
470  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
471  TYPE(alist_type), POINTER :: alist_ac, alist_bc
472  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_rrv, blocks_rrv_vrr, blocks_rv, &
473  blocks_rvr, blocks_rxrv
474  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_r_doublecom, blocks_r_rxvr, &
475  blocks_rxvr_r
476  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
477  DIMENSION(:) :: basis_set
478  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
479  TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
480 
481 !$ INTEGER(kind=omp_lock_kind), &
482 !$ ALLOCATABLE, DIMENSION(:) :: locks
483 !$ INTEGER :: lock_num, hash
484 !$ INTEGER, PARAMETER :: nlock = 501
485 
486  ppnl_present = ASSOCIATED(sap_ppnl)
487  IF (.NOT. ppnl_present) RETURN
488 
489  CALL timeset(routinen, handle)
490 
491  my_r_doublecom = .false.
492  my_r_rxvr = .false.
493  my_rxvr_r = .false.
494  my_rxrv = .false.
495  my_rrv = .false.
496  my_rv = .false.
497  my_rvr = .false.
498  my_rrv_vrr = .false.
499  IF (PRESENT(matrix_r_doublecom)) my_r_doublecom = .true.
500  IF (PRESENT(matrix_r_rxvr)) my_r_rxvr = .true.
501  IF (PRESENT(matrix_rxvr_r)) my_rxvr_r = .true.
502  IF (PRESENT(matrix_rxrv)) my_rxrv = .true.
503  IF (PRESENT(matrix_rrv)) my_rrv = .true.
504  IF (PRESENT(matrix_rv)) my_rv = .true.
505  IF (PRESENT(matrix_rvr)) my_rvr = .true.
506  IF (PRESENT(matrix_rrv_vrr)) my_rrv_vrr = .true.
507  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
508  cpabort('No dbcsr matrix provided for commutator calculation!')
509  END IF
510 
511  natom = SIZE(particle_set)
512 
513  IF (my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) THEN
514  order = 2
515  cpassert(PRESENT(ref_point)) ! need reference point for r x [r,Vnl] and [rr,Vnl]
516  ELSE IF (my_rvr .OR. my_rrv_vrr) THEN
517  order = 2
518  ELSE
519  order = 1
520  END IF
521 
522  ! When we want the double commutator [[Vnl, r], r], we also want to fix the pseudoatom
523  IF (my_r_doublecom) THEN
524  cpassert(PRESENT(pseudoatom))
525  END IF
526 
527  periodic = any(cell%perd > 0)
528  my_ref = .false.
529  IF (PRESENT(ref_point)) THEN
530  IF (.NOT. periodic) THEN
531  rf = ref_point
532  my_ref = .true.
533  ELSE ! use my_ref = False in periodic case, corresponds to distributed ref point
534  IF (order .GT. 1) THEN
535  cpwarn("Not clear how to define reference point for order > 1 in periodic cells.")
536  END IF
537  END IF
538  END IF
539 
540  nkind = SIZE(qs_kind_set)
541 
542  !sap_int needs to be shared as multiple threads need to access this
543  NULLIFY (sap_int)
544  ALLOCATE (sap_int(nkind*nkind))
545  DO i = 1, nkind*nkind
546  NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
547  sap_int(i)%nalist = 0
548  END DO
549 
550  IF (my_ref) THEN
551  ! calculate integrals <a|x^n|p>
552  CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=rf, &
553  particle_set=particle_set, cell=cell)
554  ELSE
555  CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
556  END IF
557 
558  ! *** Set up a sorting index
559  CALL sap_sort(sap_int)
560 
561  ALLOCATE (basis_set(nkind))
562  DO ikind = 1, nkind
563  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
564  IF (ASSOCIATED(orb_basis_set)) THEN
565  basis_set(ikind)%gto_basis_set => orb_basis_set
566  ELSE
567  NULLIFY (basis_set(ikind)%gto_basis_set)
568  END IF
569  END DO
570 
571  ! *** All integrals needed have been calculated and stored in sap_int
572  ! *** We now calculate the commutator matrix elements
573  CALL get_neighbor_list_set_p(neighbor_list_sets=sab_all, symmetric=do_symmetric)
574 
575 !$OMP PARALLEL &
576 !$OMP DEFAULT (NONE) &
577 !$OMP SHARED (basis_set, matrix_rv, matrix_rxrv, matrix_rrv, &
578 !$OMP matrix_rvr, matrix_rrv_vrr, matrix_r_doublecom, &
579 !$OMP sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
580 !$OMP my_rv, my_rxrv, my_rrv, my_rvr, my_rrv_vrr, &
581 !$OMP my_r_doublecom, &
582 !$OMP matrix_r_rxvr, matrix_rxvr_r, my_r_rxvr, my_rxvr_r, &
583 !$OMP pseudoatom, do_symmetric) &
584 !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
585 !$OMP iab, irow, icol, &
586 !$OMP blocks_rv, blocks_rxrv, blocks_rrv, blocks_rvr, blocks_rrv_vrr, &
587 !$OMP blocks_r_rxvr, blocks_rxvr_r, blocks_r_doublecom, &
588 !$OMP found, iac, ibc, alist_ac, alist_bc, &
589 !$OMP na, np, nb, kkind, kac, kbc, i, &
590 !$OMP go, asso_rv, asso_rxrv, asso_rrv, asso_rvr, asso_rrv_vrr, &
591 !$OMP asso_r_rxvr, asso_rxvr_r, asso_r_doublecom, hash, &
592 !$OMP acint, achint, bcint, bchint)
593 
594 !$OMP SINGLE
595 !$ ALLOCATE (locks(nlock))
596 !$OMP END SINGLE
597 
598 !$OMP DO
599 !$ DO lock_num = 1, nlock
600 !$ call omp_init_lock(locks(lock_num))
601 !$ END DO
602 !$OMP END DO
603 
604 !$OMP DO SCHEDULE(GUIDED)
605 
606  DO slot = 1, sab_all(1)%nl_size
607 
608  ikind = sab_all(1)%nlist_task(slot)%ikind
609  jkind = sab_all(1)%nlist_task(slot)%jkind
610  iatom = sab_all(1)%nlist_task(slot)%iatom
611  jatom = sab_all(1)%nlist_task(slot)%jatom
612  cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
613  rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
614 
615  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
616  IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
617  iab = ikind + nkind*(jkind - 1)
618 
619  IF (do_symmetric) THEN
620  IF (iatom <= jatom) THEN
621  irow = iatom
622  icol = jatom
623  ELSE
624  irow = jatom
625  icol = iatom
626  END IF
627  ELSE
628  irow = iatom
629  icol = jatom
630  END IF
631 
632  ! allocate blocks
633  IF (my_rv) THEN
634  ALLOCATE (blocks_rv(3))
635  END IF
636  IF (my_rxrv) THEN
637  ALLOCATE (blocks_rxrv(3))
638  END IF
639  IF (my_rrv) THEN
640  ALLOCATE (blocks_rrv(6))
641  END IF
642  IF (my_rvr) THEN
643  ALLOCATE (blocks_rvr(6))
644  END IF
645  IF (my_rrv_vrr) THEN
646  ALLOCATE (blocks_rrv_vrr(6))
647  END IF
648  IF (my_r_rxvr) THEN
649  ALLOCATE (blocks_r_rxvr(3, 3))
650  END IF
651 
652  IF (my_rxvr_r) THEN
653  ALLOCATE (blocks_rxvr_r(3, 3))
654  END IF
655 
656  IF (my_r_doublecom) THEN
657  ALLOCATE (blocks_r_doublecom(3, 3))
658  END IF
659 
660  ! get blocks
661  IF (my_rv) THEN
662  DO ind = 1, 3
663  CALL dbcsr_get_block_p(matrix_rv(ind)%matrix, irow, icol, blocks_rv(ind)%block, found)
664  END DO
665  END IF
666 
667  IF (my_rxrv) THEN
668  DO ind = 1, 3
669  CALL dbcsr_get_block_p(matrix_rxrv(ind)%matrix, irow, icol, blocks_rxrv(ind)%block, found)
670  blocks_rxrv(ind)%block(:, :) = 0._dp
671  END DO
672  END IF
673 
674  IF (my_rrv) THEN
675  DO ind = 1, 6
676  CALL dbcsr_get_block_p(matrix_rrv(ind)%matrix, irow, icol, blocks_rrv(ind)%block, found)
677  END DO
678  END IF
679 
680  IF (my_rvr) THEN
681  DO ind = 1, 6
682  CALL dbcsr_get_block_p(matrix_rvr(ind)%matrix, irow, icol, blocks_rvr(ind)%block, found)
683  END DO
684  END IF
685 
686  IF (my_rrv_vrr) THEN
687  DO ind = 1, 6
688  CALL dbcsr_get_block_p(matrix_rrv_vrr(ind)%matrix, irow, icol, blocks_rrv_vrr(ind)%block, found)
689  END DO
690  END IF
691 
692  IF (my_r_rxvr) THEN
693  DO ind = 1, 3
694  DO ind2 = 1, 3
695  CALL dbcsr_get_block_p(matrix_r_rxvr(ind, ind2)%matrix, irow, icol, &
696  blocks_r_rxvr(ind, ind2)%block, found)
697  blocks_r_rxvr(ind, ind2)%block(:, :) = 0._dp
698  END DO
699  END DO
700  END IF
701 
702  IF (my_rxvr_r) THEN
703  DO ind = 1, 3
704  DO ind2 = 1, 3
705  CALL dbcsr_get_block_p(matrix_rxvr_r(ind, ind2)%matrix, irow, icol, &
706  blocks_rxvr_r(ind, ind2)%block, found)
707  blocks_rxvr_r(ind, ind2)%block(:, :) = 0._dp
708  END DO
709  END DO
710  END IF
711 
712  IF (my_r_doublecom) THEN
713  DO ind = 1, 3
714  DO ind2 = 1, 3
715  CALL dbcsr_get_block_p(matrix_r_doublecom(ind, ind2)%matrix, irow, icol, &
716  blocks_r_doublecom(ind, ind2)%block, found)
717  blocks_r_doublecom(ind, ind2)%block(:, :) = 0._dp
718  END DO
719  END DO
720  END IF
721 
722  ! check whether all blocks are associated
723  go = .true.
724  IF (my_rv) THEN
725  asso_rv = (ASSOCIATED(blocks_rv(1)%block) .AND. ASSOCIATED(blocks_rv(2)%block) .AND. &
726  ASSOCIATED(blocks_rv(3)%block))
727  go = go .AND. asso_rv
728  END IF
729 
730  IF (my_rxrv) THEN
731  asso_rxrv = (ASSOCIATED(blocks_rxrv(1)%block) .AND. ASSOCIATED(blocks_rxrv(2)%block) .AND. &
732  ASSOCIATED(blocks_rxrv(3)%block))
733  go = go .AND. asso_rxrv
734  END IF
735 
736  IF (my_rrv) THEN
737  asso_rrv = (ASSOCIATED(blocks_rrv(1)%block) .AND. ASSOCIATED(blocks_rrv(2)%block) .AND. &
738  ASSOCIATED(blocks_rrv(3)%block) .AND. ASSOCIATED(blocks_rrv(4)%block) .AND. &
739  ASSOCIATED(blocks_rrv(5)%block) .AND. ASSOCIATED(blocks_rrv(6)%block))
740  go = go .AND. asso_rrv
741  END IF
742 
743  IF (my_rvr) THEN
744  asso_rvr = (ASSOCIATED(blocks_rvr(1)%block) .AND. ASSOCIATED(blocks_rvr(2)%block) .AND. &
745  ASSOCIATED(blocks_rvr(3)%block) .AND. ASSOCIATED(blocks_rvr(4)%block) .AND. &
746  ASSOCIATED(blocks_rvr(5)%block) .AND. ASSOCIATED(blocks_rvr(6)%block))
747  go = go .AND. asso_rvr
748  END IF
749 
750  IF (my_rrv_vrr) THEN
751  asso_rrv_vrr = (ASSOCIATED(blocks_rrv_vrr(1)%block) .AND. ASSOCIATED(blocks_rrv_vrr(2)%block) .AND. &
752  ASSOCIATED(blocks_rrv_vrr(3)%block) .AND. ASSOCIATED(blocks_rrv_vrr(4)%block) .AND. &
753  ASSOCIATED(blocks_rrv_vrr(5)%block) .AND. ASSOCIATED(blocks_rrv_vrr(6)%block))
754  go = go .AND. asso_rrv_vrr
755  END IF
756 
757  IF (my_r_rxvr) THEN
758  asso_r_rxvr = .true.
759  DO ind = 1, 3
760  DO ind2 = 1, 3
761  asso_r_rxvr = asso_r_rxvr .AND. ASSOCIATED(blocks_r_rxvr(ind, ind2)%block)
762  END DO
763  END DO
764  go = go .AND. asso_r_rxvr
765  END IF
766 
767  IF (my_rxvr_r) THEN
768  asso_rxvr_r = .true.
769  DO ind = 1, 3
770  DO ind2 = 1, 3
771  asso_rxvr_r = asso_rxvr_r .AND. ASSOCIATED(blocks_rxvr_r(ind, ind2)%block)
772  END DO
773  END DO
774  go = go .AND. asso_rxvr_r
775  END IF
776 
777  IF (my_r_doublecom) THEN
778  asso_r_doublecom = .true.
779  DO ind = 1, 3
780  DO ind2 = 1, 3
781  asso_r_doublecom = asso_r_doublecom .AND. ASSOCIATED(blocks_r_doublecom(ind, ind2)%block)
782  END DO
783  END DO
784  go = go .AND. asso_r_doublecom
785  END IF
786 
787  ! loop over all kinds for projector atom
788  ! < iatom | katom > h < katom | jatom >
789  IF (go) THEN
790  DO kkind = 1, nkind
791  iac = ikind + nkind*(kkind - 1)
792  ibc = jkind + nkind*(kkind - 1)
793  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
794  IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
795  CALL get_alist(sap_int(iac), alist_ac, iatom)
796  CALL get_alist(sap_int(ibc), alist_bc, jatom)
797  IF (.NOT. ASSOCIATED(alist_ac)) cycle
798  IF (.NOT. ASSOCIATED(alist_bc)) cycle
799  DO kac = 1, alist_ac%nclist
800  DO kbc = 1, alist_bc%nclist
801  IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
802  IF (PRESENT(pseudoatom)) THEN
803  IF (alist_ac%clist(kac)%catom /= pseudoatom) cycle
804  END IF
805 
806  IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
807  IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
808  acint => alist_ac%clist(kac)%acint
809  bcint => alist_bc%clist(kbc)%acint
810  achint => alist_ac%clist(kac)%achint
811  bchint => alist_bc%clist(kbc)%achint
812  na = SIZE(acint, 1)
813  np = SIZE(acint, 2)
814  nb = SIZE(bcint, 1)
815 !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
816 !$ CALL omp_set_lock(locks(hash))
817  IF (my_rv) THEN
818  ! r*Vnl
819  ! with LAPACK
820  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 2), na, &
821  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(1)%block, SIZE(blocks_rv(1)%block, 1)) ! xV
822  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 3), na, &
823  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(2)%block, SIZE(blocks_rv(2)%block, 1)) ! yV
824  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 4), na, &
825  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(3)%block, SIZE(blocks_rv(3)%block, 1)) ! zV
826  IF (iatom <= jatom) THEN
827  ! with MATMUL
828  blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
829  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1))) ! xV
830  blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) + &
831  matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1))) ! yV
832  blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) + &
833  matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1))) ! zV
834  ELSE
835  blocks_rv(1)%block(1:nb, 1:na) = blocks_rv(1)%block(1:nb, 1:na) + &
836  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))
837  blocks_rv(2)%block(1:nb, 1:na) = blocks_rv(2)%block(1:nb, 1:na) + &
838  matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))
839  blocks_rv(3)%block(1:nb, 1:na) = blocks_rv(3)%block(1:nb, 1:na) + &
840  matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))
841  END IF
842  ! -Vnl r
843  ! with LAPACK
844  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
845  ! bcint(1, 1, 2), nb, 1.0_dp, blocks_rv(1)%block, SIZE(blocks_rv(1)%block, 1)) ! -Vx
846  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
847  ! bcint(1, 1, 3), nb, 1.0_dp, blocks_rv(2)%block, SIZE(blocks_rv(2)%block, 1)) ! -Vy
848  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
849  ! bcint(1, 1, 4), nb, 1.0_dp, blocks_rv(3)%block, SIZE(blocks_rv(3)%block, 1)) ! -Vz
850  ! with MATMUL
851  IF (iatom <= jatom) THEN
852  blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) - &
853  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) ! -Vx
854  blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) - &
855  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) ! -Vy
856  blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) - &
857  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) ! -Vz
858  ELSE
859  blocks_rv(1)%block(1:nb, 1:na) = blocks_rv(1)%block(1:nb, 1:na) - &
860  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2)))
861  blocks_rv(2)%block(1:nb, 1:na) = blocks_rv(2)%block(1:nb, 1:na) - &
862  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3)))
863  blocks_rv(3)%block(1:nb, 1:na) = blocks_rv(3)%block(1:nb, 1:na) - &
864  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4)))
865  END IF
866 
867  END IF
868 
869  IF (my_rxrv) THEN
870  ! x-component (y [z,Vnl] - z [y, Vnl])
871  ! with LAPACK
872  ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 9), na, &
873  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! yzV
874  ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 3), na, &
875  ! bcint(1, 1, 4), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! -yVz
876  ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 9), na, &
877  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! -zyV
878  ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 4), na, &
879  ! bcint(1, 1, 3), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! zVy
880  ! with MATMUL
881  IF (iatom <= jatom) THEN
882  blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
883  matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1))) ! yzV
884  blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
885  matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4))) ! -yVz
886  blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
887  matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1))) ! -zyV
888  blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
889  matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 3))) ! zVy
890  ELSE
891  blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
892  matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1))) ! yzV
893  blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
894  matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4))) ! -yVz
895  blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
896  matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1))) ! -zyV
897  blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
898  matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 3))) ! zVy
899  END IF
900 
901  ! y-component (z [x,Vnl] - x [z, Vnl])
902  ! with LAPACK
903  ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 7), na, &
904  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! zxV
905  ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 4), na, &
906  ! bcint(1, 1, 2), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! -zVx
907  ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 7), na, &
908  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! -xzV
909  ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 2), na, &
910  ! bcint(1, 1, 4), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! xVz
911  ! with MATMUL
912  IF (iatom <= jatom) THEN
913  blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
914  matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1))) ! zxV
915  blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
916  matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 2))) ! -zVx
917  blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
918  matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1))) ! -xzV
919  blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
920  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4))) ! xVz
921  ELSE
922  blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
923  matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1))) ! zxV
924  blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
925  matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 2))) ! -zVx
926  blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
927  matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1))) ! -xzV
928  blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
929  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4))) ! xVz
930  END IF
931 
932  ! z-component (x [y,Vnl] - y [x, Vnl])
933  ! with LAPACK
934  ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 6), na, &
935  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! xyV
936  ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 2), na, &
937  ! bcint(1, 1, 3), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! -xVy
938  ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 6), na, &
939  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! -yxV
940  ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 3), na, &
941  ! bcint(1, 1, 2), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! yVx
942  ! with MATMUL
943  IF (iatom <= jatom) THEN
944  blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
945  matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1))) ! xyV
946  blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
947  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3))) ! -xVy
948  blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
949  matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1))) ! -yxV
950  blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
951  matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 2))) ! zVx
952  ELSE
953  blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
954  matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1))) ! xyV
955  blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
956  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3))) ! -xVy
957  blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
958  matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1))) ! -yxV
959  blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
960  matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 2))) ! zVx
961  END IF
962  END IF
963 
964  IF (my_rrv) THEN
965  ! r_alpha * r_beta * Vnl
966  ! with LAPACK
967  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 5), na, &
968  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(1)%block, SIZE(blocks_rrv(1)%block, 1)) ! xxV
969  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 6), na, &
970  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(2)%block, SIZE(blocks_rrv(2)%block, 1)) ! xyV
971  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 7), na, &
972  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(3)%block, SIZE(blocks_rrv(3)%block, 1)) ! xzV
973  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 8), na, &
974  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(4)%block, SIZE(blocks_rrv(4)%block, 1)) ! yyV
975  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 9), na, &
976  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(5)%block, SIZE(blocks_rrv(5)%block, 1)) ! yzV
977  ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 10), na, &
978  ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(6)%block, SIZE(blocks_rrv(6)%block, 1)) ! zzV
979  ! with MATMUL
980  IF (iatom <= jatom) THEN
981  blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) + &
982  matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1))) ! xxV
983  blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) + &
984  matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1))) ! xyV
985  blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) + &
986  matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1))) ! xzV
987  blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) + &
988  matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1))) ! yyV
989  blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) + &
990  matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1))) ! yzV
991  blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) + &
992  matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1))) ! zzV
993  ELSE
994  blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) + &
995  matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1))) ! xxV
996  blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) + &
997  matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1))) ! xyV
998  blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) + &
999  matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1))) ! xzV
1000  blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) + &
1001  matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1))) ! yyV
1002  blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) + &
1003  matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1))) ! yzV
1004  blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) + &
1005  matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1))) ! zzV
1006  END IF
1007 
1008  ! - Vnl * r_alpha * r_beta
1009  ! with LAPACK
1010  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1011  ! bcint(1, 1, 5), nb, 1.0_dp, blocks_rrv(1)%block, SIZE(blocks_rrv(1)%block, 1)) ! Vxx
1012  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1013  ! bcint(1, 1, 6), nb, 1.0_dp, blocks_rrv(2)%block, SIZE(blocks_rrv(2)%block, 1)) ! Vxy
1014  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1015  ! bcint(1, 1, 7), nb, 1.0_dp, blocks_rrv(3)%block, SIZE(blocks_rrv(3)%block, 1)) ! Vxz
1016  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1017  ! bcint(1, 1, 8), nb, 1.0_dp, blocks_rrv(4)%block, SIZE(blocks_rrv(4)%block, 1)) ! Vyy
1018  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1019  ! bcint(1, 1, 9), nb, 1.0_dp, blocks_rrv(5)%block, SIZE(blocks_rrv(5)%block, 1)) ! Vyz
1020  ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1021  ! bcint(1, 1, 10), nb, 1.0_dp, blocks_rrv(6)%block, SIZE(blocks_rrv(6)%block, 1)) ! Vzz
1022  ! with MATMUL
1023  IF (iatom <= jatom) THEN
1024  blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) - &
1025  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5))) ! -Vxx
1026  blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) - &
1027  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6))) ! -Vxy
1028  blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) - &
1029  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7))) ! -Vxz
1030  blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) - &
1031  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8))) ! -Vyy
1032  blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) - &
1033  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9))) ! -Vyz
1034  blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) - &
1035  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10))) ! -Vzz
1036  ELSE
1037  blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) - &
1038  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5))) ! -Vxx
1039  blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) - &
1040  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6))) ! -Vxy
1041  blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) - &
1042  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7))) ! -Vxz
1043  blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) - &
1044  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8))) ! -Vyy
1045  blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) - &
1046  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9))) ! -Vyz
1047  blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) - &
1048  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10))) ! -Vzz
1049  END IF
1050  END IF
1051 
1052  IF (my_rvr) THEN
1053  ! r_alpha * Vnl * r_beta
1054  IF (iatom <= jatom) THEN
1055  blocks_rvr(1)%block(1:na, 1:nb) = blocks_rvr(1)%block(1:na, 1:nb) + &
1056  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 2))) ! xVx
1057  blocks_rvr(2)%block(1:na, 1:nb) = blocks_rvr(2)%block(1:na, 1:nb) + &
1058  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3))) ! xVy
1059  blocks_rvr(3)%block(1:na, 1:nb) = blocks_rvr(3)%block(1:na, 1:nb) + &
1060  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4))) ! xVz
1061  blocks_rvr(4)%block(1:na, 1:nb) = blocks_rvr(4)%block(1:na, 1:nb) + &
1062  matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 3))) ! yVy
1063  blocks_rvr(5)%block(1:na, 1:nb) = blocks_rvr(5)%block(1:na, 1:nb) + &
1064  matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4))) ! yVz
1065  blocks_rvr(6)%block(1:na, 1:nb) = blocks_rvr(6)%block(1:na, 1:nb) + &
1066  matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 4))) ! zVz
1067  ELSE
1068  blocks_rvr(1)%block(1:nb, 1:na) = blocks_rvr(1)%block(1:nb, 1:na) + &
1069  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 2))) ! xVx
1070  blocks_rvr(2)%block(1:nb, 1:na) = blocks_rvr(2)%block(1:nb, 1:na) + &
1071  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3))) ! xVy
1072  blocks_rvr(3)%block(1:nb, 1:na) = blocks_rvr(3)%block(1:nb, 1:na) + &
1073  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4))) ! xVz
1074  blocks_rvr(4)%block(1:nb, 1:na) = blocks_rvr(4)%block(1:nb, 1:na) + &
1075  matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 3))) ! yVy
1076  blocks_rvr(5)%block(1:nb, 1:na) = blocks_rvr(5)%block(1:nb, 1:na) + &
1077  matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4))) ! yVz
1078  blocks_rvr(6)%block(1:nb, 1:na) = blocks_rvr(6)%block(1:nb, 1:na) + &
1079  matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 4))) ! zVz
1080  END IF
1081  END IF
1082 
1083  IF (my_rrv_vrr) THEN
1084  ! r_alpha * r_beta * Vnl
1085  IF (iatom <= jatom) THEN
1086  blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
1087  matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1))) ! xxV
1088  blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
1089  matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1))) ! xyV
1090  blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
1091  matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1))) ! xzV
1092  blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
1093  matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1))) ! yyV
1094  blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
1095  matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1))) ! yzV
1096  blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
1097  matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1))) ! zzV
1098  ELSE
1099  blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
1100  matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1))) ! xxV
1101  blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
1102  matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1))) ! xyV
1103  blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
1104  matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1))) ! xzV
1105  blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
1106  matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1))) ! yyV
1107  blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
1108  matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1))) ! yzV
1109  blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
1110  matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1))) ! zzV
1111  END IF
1112  ! + Vnl * r_alpha * r_beta
1113  IF (iatom <= jatom) THEN
1114  blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
1115  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5))) ! +Vxx
1116  blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
1117  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6))) ! +Vxy
1118  blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
1119  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7))) ! +Vxz
1120  blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
1121  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8))) ! +Vyy
1122  blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
1123  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9))) ! +Vyz
1124  blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
1125  matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10))) ! +Vzz
1126  ELSE
1127  blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
1128  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5))) ! +Vxx
1129  blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
1130  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6))) ! +Vxy
1131  blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
1132  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7))) ! +Vxz
1133  blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
1134  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8))) ! +Vyy
1135  blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
1136  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9))) ! +Vyz
1137  blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
1138  matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10))) ! +Vzz
1139  END IF
1140  END IF
1141 
1142  ! The indices are stored in i_1, i_x, ..., i_zzz
1143  ! matrix_r_rxvr(alpha, beta)
1144  ! = sum_(gamma delta) epsilon_(alpha gamma delta)
1145  ! (r_beta * r_gamma * V_nl * r_delta - r_beta * r_gamma * r_delta * V_nl)
1146  ! = sum_(gamma delta) epsilon_(alpha gamma delta) r_beta * r_gamma * V_nl * r_delta
1147 
1148  ! TODO: is this set to zero before?
1149  IF (my_r_rxvr) THEN
1150  ! beta = 1
1151  ! matrix_r_rxvr(x, x) = x * y * V_nl * z - x * z * V_nl * y
1152  blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
1153  blocks_r_rxvr(1, 1)%block(1:na, 1:nb) + &
1154  matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
1155  blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
1156  blocks_r_rxvr(1, 1)%block(1:na, 1:nb) - &
1157  matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
1158 
1159  ! matrix_r_rxvr(y, x) = x * z * V_nl * x - x * x * V_nl * z
1160  blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
1161  blocks_r_rxvr(2, 1)%block(1:na, 1:nb) + &
1162  matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_x)))
1163  blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
1164  blocks_r_rxvr(2, 1)%block(1:na, 1:nb) - &
1165  matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
1166 
1167  ! matrix_r_rxvr(z, x) = x * x * V_nl * y - x * y * V_nl * x
1168  blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
1169  blocks_r_rxvr(3, 1)%block(1:na, 1:nb) + &
1170  matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
1171  blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
1172  blocks_r_rxvr(3, 1)%block(1:na, 1:nb) - &
1173  matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_x)))
1174 
1175  ! beta = 2
1176  ! matrix_r_rxvr(x, y) = y * y * V_nl * z - y * z * V_nl * y
1177  blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
1178  blocks_r_rxvr(1, 2)%block(1:na, 1:nb) + &
1179  matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
1180  blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
1181  blocks_r_rxvr(1, 2)%block(1:na, 1:nb) - &
1182  matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_y)))
1183 
1184  ! matrix_r_rxvr(y, y) = y * z * V_nl * x - y * x * V_nl * z
1185  blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
1186  blocks_r_rxvr(2, 2)%block(1:na, 1:nb) + &
1187  matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
1188  blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
1189  blocks_r_rxvr(2, 2)%block(1:na, 1:nb) - &
1190  matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
1191 
1192  ! matrix_r_rxvr(z, y) = y * x * V_nl * y - y * y * V_nl * x
1193  blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
1194  blocks_r_rxvr(3, 2)%block(1:na, 1:nb) + &
1195  matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_y)))
1196  blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
1197  blocks_r_rxvr(3, 2)%block(1:na, 1:nb) - &
1198  matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
1199 
1200  ! beta = 3
1201  ! matrix_r_rxvr(x, z) = z * y * V_nl * z - z * z * V_nl * y
1202  blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
1203  blocks_r_rxvr(1, 3)%block(1:na, 1:nb) + &
1204  matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_z)))
1205  blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
1206  blocks_r_rxvr(1, 3)%block(1:na, 1:nb) - &
1207  matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
1208 
1209  ! matrix_r_rxvr(y, z) = z * z * V_nl * x - z * x * V_nl * z
1210  blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
1211  blocks_r_rxvr(2, 3)%block(1:na, 1:nb) + &
1212  matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
1213  blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
1214  blocks_r_rxvr(2, 3)%block(1:na, 1:nb) - &
1215  matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_z)))
1216 
1217  ! matrix_r_rxvr(z, z) = z * x * V_nl * y - z * y * V_nl * x
1218  blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
1219  blocks_r_rxvr(3, 3)%block(1:na, 1:nb) + &
1220  matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
1221  blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
1222  blocks_r_rxvr(3, 3)%block(1:na, 1:nb) - &
1223  matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
1224 
1225  END IF ! my_r_rxvr
1226 
1227  ! The indices are stored in i_1, i_x, ..., i_zzz
1228  ! This will put into blocks_rxvr_r
1229  ! matrix_rxvr_r(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
1230  ! r_gamma * V_nl * r_delta * r_beta
1231  IF (my_rxvr_r) THEN
1232  ! beta = 1
1233  ! matrix_rxvr_r(x, x) = yV zx - zV yx
1234  blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
1235  blocks_rxvr_r(1, 1)%block(1:na, 1:nb) + &
1236  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
1237  blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
1238  blocks_rxvr_r(1, 1)%block(1:na, 1:nb) - &
1239  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
1240 
1241  ! matrix_rxvr_r(y, x) = zV xx - xV zx
1242  blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
1243  blocks_rxvr_r(2, 1)%block(1:na, 1:nb) + &
1244  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
1245  blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
1246  blocks_rxvr_r(2, 1)%block(1:na, 1:nb) - &
1247  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zx)))
1248 
1249  ! matrix_rxvr_r(z, x) = xV yx - yV xx
1250  blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
1251  blocks_rxvr_r(3, 1)%block(1:na, 1:nb) + &
1252  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yx)))
1253  blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
1254  blocks_rxvr_r(3, 1)%block(1:na, 1:nb) - &
1255  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
1256 
1257  ! beta = 2
1258  ! matrix_rxvr_r(x, y) = yV zy - zV yy
1259  blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
1260  blocks_rxvr_r(1, 2)%block(1:na, 1:nb) + &
1261  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zy)))
1262  blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
1263  blocks_rxvr_r(1, 2)%block(1:na, 1:nb) - &
1264  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
1265 
1266  ! matrix_rxvr_r(y, y) = zV xy - xV zy
1267  blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
1268  blocks_rxvr_r(2, 2)%block(1:na, 1:nb) + &
1269  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
1270  blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
1271  blocks_rxvr_r(2, 2)%block(1:na, 1:nb) - &
1272  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
1273 
1274  ! matrix_rxvr_r(z, y) = xV yy - yV xy
1275  blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
1276  blocks_rxvr_r(3, 2)%block(1:na, 1:nb) + &
1277  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
1278  blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
1279  blocks_rxvr_r(3, 2)%block(1:na, 1:nb) - &
1280  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xy)))
1281 
1282  ! beta = 3
1283  ! matrix_rxvr_r(x, z) = yV zz - zV yz
1284  blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
1285  blocks_rxvr_r(1, 3)%block(1:na, 1:nb) + &
1286  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
1287  blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
1288  blocks_rxvr_r(1, 3)%block(1:na, 1:nb) - &
1289  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yz)))
1290 
1291  ! matrix_rxvr_r(y, z) = zV xz - xV zz
1292  blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
1293  blocks_rxvr_r(2, 3)%block(1:na, 1:nb) + &
1294  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xz)))
1295  blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
1296  blocks_rxvr_r(2, 3)%block(1:na, 1:nb) - &
1297  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
1298 
1299  ! matrix_rxvr_r(z, z) = xV yz - yV xz
1300  blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
1301  blocks_rxvr_r(3, 3)%block(1:na, 1:nb) + &
1302  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
1303  blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
1304  blocks_rxvr_r(3, 3)%block(1:na, 1:nb) - &
1305  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
1306 
1307  END IF ! my_rxvr_r
1308 
1309  ! matrix_r_doublecom(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
1310  ! gamma V^pseudoatom beta delta - gamma beta V^pseudoatom delta
1311 
1312  IF (my_r_doublecom) THEN
1313  ! beta = 1
1314  ! matrix_r_doublecom(x, x) = yV xz - zV xy - yxV z + zxV y
1315  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1316  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
1317  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
1318  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1319  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
1320  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
1321  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1322  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
1323  matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
1324  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1325  blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
1326  matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
1327 
1328  ! matrix_r_doublecom(y, x) = zV xx - xV xz - zxV x + xxV z
1329  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1330  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
1331  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
1332  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1333  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
1334  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xz)))
1335  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1336  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
1337  matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_x)))
1338  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1339  blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
1340  matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
1341 
1342  ! matrix_r_doublecom(z, x) = xV xy - yV xx - xxV y + yxV x
1343  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1344  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1345  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xy)))
1346  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1347  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1348  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
1349  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1350  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1351  matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
1352  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1353  blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1354  matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_x)))
1355 
1356  ! beta = 2
1357  ! matrix_r_doublecom(x, y) = yV yz - zV yy - yyV z + zyV y
1358  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1359  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1360  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yz)))
1361  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1362  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1363  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
1364  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1365  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1366  matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
1367  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1368  blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1369  matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_y)))
1370 
1371  ! matrix_r_doublecom(y, y) = zV yx - xV yz - zyV x + xyV z
1372  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1373  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1374  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
1375  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1376  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1377  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
1378  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1379  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1380  matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
1381  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1382  blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1383  matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
1384 
1385  ! matrix_r_doublecom(z, y) = xV yy - yV yx - xyV y + yyV x
1386  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1387  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1388  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
1389  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1390  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1391  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yx)))
1392  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1393  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1394  matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_y)))
1395  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1396  blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1397  matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
1398 
1399  ! beta = 3
1400  ! matrix_r_doublecom(x, z) = yV zz - zV zy - yzV z + zzV y
1401  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1402  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1403  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
1404  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1405  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1406  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zy)))
1407  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1408  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1409  matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_z)))
1410  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1411  blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1412  matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
1413 
1414  ! matrix_r_doublecom(y, z) = zV zx - xV zz - zzV x + xzV z
1415  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1416  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1417  matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zx)))
1418  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1419  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1420  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
1421  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1422  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1423  matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
1424  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1425  blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1426  matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_z)))
1427 
1428  ! matrix_r_doublecom(z, z) = xV zy - yV zx - xzV y + yzV x
1429  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1430  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1431  matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
1432  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1433  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1434  matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
1435  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1436  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1437  matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
1438  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1439  blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1440  matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
1441 
1442  END IF ! my_r_doublecom
1443 !$ CALL omp_unset_lock(locks(hash))
1444  EXIT ! We have found a match and there can be only one single match
1445  END IF
1446  END DO
1447  END DO
1448  END DO
1449  END IF
1450  IF (my_rv) THEN
1451  DO ind = 1, 3
1452  NULLIFY (blocks_rv(ind)%block)
1453  END DO
1454  DEALLOCATE (blocks_rv)
1455  END IF
1456  IF (my_rxrv) THEN
1457  DO ind = 1, 3
1458  NULLIFY (blocks_rxrv(ind)%block)
1459  END DO
1460  DEALLOCATE (blocks_rxrv)
1461  END IF
1462  IF (my_rrv) THEN
1463  DO ind = 1, 6
1464  NULLIFY (blocks_rrv(ind)%block)
1465  END DO
1466  DEALLOCATE (blocks_rrv)
1467  END IF
1468  IF (my_rvr) THEN
1469  DO ind = 1, 6
1470  NULLIFY (blocks_rvr(ind)%block)
1471  END DO
1472  DEALLOCATE (blocks_rvr)
1473  END IF
1474  IF (my_rrv_vrr) THEN
1475  DO ind = 1, 6
1476  NULLIFY (blocks_rrv_vrr(ind)%block)
1477  END DO
1478  DEALLOCATE (blocks_rrv_vrr)
1479  END IF
1480  IF (my_r_rxvr) THEN
1481  DO ind = 1, 3
1482  DO ind2 = 1, 3
1483  NULLIFY (blocks_r_rxvr(ind, ind2)%block)
1484  END DO
1485  END DO
1486  DEALLOCATE (blocks_r_rxvr)
1487  END IF
1488  IF (my_rxvr_r) THEN
1489  DO ind = 1, 3
1490  DO ind2 = 1, 3
1491  NULLIFY (blocks_rxvr_r(ind, ind2)%block)
1492  END DO
1493  END DO
1494  DEALLOCATE (blocks_rxvr_r)
1495  END IF
1496  IF (my_r_doublecom) THEN
1497  DO ind = 1, 3
1498  DO ind2 = 1, 3
1499  NULLIFY (blocks_r_doublecom(ind, ind2)%block)
1500  END DO
1501  END DO
1502  DEALLOCATE (blocks_r_doublecom)
1503  END IF
1504  END DO
1505 
1506 !$OMP DO
1507 !$ DO lock_num = 1, nlock
1508 !$ call omp_destroy_lock(locks(lock_num))
1509 !$ END DO
1510 !$OMP END DO
1511 
1512 !$OMP SINGLE
1513 !$ DEALLOCATE (locks)
1514 !$OMP END SINGLE NOWAIT
1515 
1516 !$OMP END PARALLEL
1517 
1518  CALL release_sap_int(sap_int)
1519 
1520  DEALLOCATE (basis_set)
1521 
1522  CALL timestop(handle)
1523 
1524  END SUBROUTINE build_com_mom_nl
1525 
1526 ! **************************************************************************************************
1527 !> \brief calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
1528 !> \param qs_kind_set ...
1529 !> \param sab_all ...
1530 !> \param sap_ppnl ...
1531 !> \param eps_ppnl ...
1532 !> \param particle_set ...
1533 !> \param matrix_mag_nl ...
1534 !> \param refpoint ...
1535 !> \param cell ...
1536 ! **************************************************************************************************
1537  SUBROUTINE build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
1538 
1539  TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1540  POINTER :: qs_kind_set
1541  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1542  INTENT(IN), POINTER :: sab_all, sap_ppnl
1543  REAL(kind=dp), INTENT(IN) :: eps_ppnl
1544  TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1545  POINTER :: particle_set
1546  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1547  POINTER :: matrix_mag_nl
1548  REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
1549  TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1550 
1551  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_nl_mag'
1552 
1553  INTEGER :: handle, iab, iac, iatom, ibc, icol, &
1554  ikind, ind, irow, jatom, jkind, kac, &
1555  kbc, kkind, na, natom, nb, nkind, np, &
1556  order, slot
1557  INTEGER, DIMENSION(3) :: cell_b
1558  LOGICAL :: found, go, my_ref, ppnl_present
1559  REAL(kind=dp), DIMENSION(3) :: r_b, r_ps, rab
1560  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1561  TYPE(alist_type), POINTER :: alist_ac, alist_bc
1562  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1563  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1564  DIMENSION(:) :: basis_set
1565  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1566  TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1567 
1568 !$ INTEGER(kind=omp_lock_kind), &
1569 !$ ALLOCATABLE, DIMENSION(:) :: locks
1570 !$ INTEGER :: lock_num, hash
1571 !$ INTEGER, PARAMETER :: nlock = 501
1572 
1573  ppnl_present = ASSOCIATED(sap_ppnl)
1574  IF (.NOT. ppnl_present) RETURN
1575 
1576  CALL timeset(routinen, handle)
1577 
1578  my_ref = .false.
1579  IF (PRESENT(refpoint)) THEN
1580  my_ref = .true.
1581  cpassert(PRESENT(cell))
1582  END IF
1583 
1584  natom = SIZE(particle_set)
1585  nkind = SIZE(qs_kind_set)
1586 
1587  ! allocate integral storage
1588  NULLIFY (sap_int)
1589  ALLOCATE (sap_int(nkind*nkind))
1590  DO ind = 1, nkind*nkind
1591  NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1592  sap_int(ind)%nalist = 0
1593  END DO
1594 
1595  ! build integrals over GTO + projector functions, refpoint actually
1596  order = 1 ! only need first moments (x, y, z)
1597  ! refpoint actually does not matter in this case, i. e. (order = 1 .and. commutator)
1598  IF (my_ref) THEN
1599  CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=refpoint, &
1600  particle_set=particle_set, cell=cell)
1601  ELSE
1602  CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1603  END IF
1604 
1605  CALL sap_sort(sap_int)
1606 
1607  ! get access to basis sets
1608  ALLOCATE (basis_set(nkind))
1609  DO ikind = 1, nkind
1610  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1611  IF (ASSOCIATED(orb_basis_set)) THEN
1612  basis_set(ikind)%gto_basis_set => orb_basis_set
1613  ELSE
1614  NULLIFY (basis_set(ikind)%gto_basis_set)
1615  END IF
1616  END DO
1617 
1618 !$OMP PARALLEL &
1619 !$OMP DEFAULT (NONE) &
1620 !$OMP SHARED (basis_set, matrix_mag_nl, sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
1621 !$OMP particle_set, my_ref, refpoint) &
1622 !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1623 !$OMP iab, irow, icol, blocks_mag, r_ps, r_b, go, hash, &
1624 !$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
1625 !$OMP achint, bchint, na, np, nb, kkind, kac, kbc)
1626 
1627 !$OMP SINGLE
1628 !$ ALLOCATE (locks(nlock))
1629 !$OMP END SINGLE
1630 
1631 !$OMP DO
1632 !$ DO lock_num = 1, nlock
1633 !$ call omp_init_lock(locks(lock_num))
1634 !$ END DO
1635 !$OMP END DO
1636 
1637 !$OMP DO SCHEDULE(GUIDED)
1638  DO slot = 1, sab_all(1)%nl_size
1639  ! get indices
1640  ikind = sab_all(1)%nlist_task(slot)%ikind
1641  jkind = sab_all(1)%nlist_task(slot)%jkind
1642  iatom = sab_all(1)%nlist_task(slot)%iatom
1643  jatom = sab_all(1)%nlist_task(slot)%jatom
1644  cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1645  rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1646 
1647  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1648  IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1649  iab = ikind + nkind*(jkind - 1)
1650 
1651  IF (iatom <= jatom) THEN
1652  irow = iatom
1653  icol = jatom
1654  ELSE
1655  irow = jatom
1656  icol = iatom
1657  END IF
1658 
1659  ! get blocks
1660  ALLOCATE (blocks_mag(3))
1661  DO ind = 1, 3
1662  CALL dbcsr_get_block_p(matrix_mag_nl(ind)%matrix, irow, icol, blocks_mag(ind)%block, found)
1663  END DO
1664 
1665  go = (ASSOCIATED(blocks_mag(1)%block) .AND. ASSOCIATED(blocks_mag(2)%block) .AND. ASSOCIATED(blocks_mag(3)%block))
1666 
1667  IF (go) THEN
1668  DO kkind = 1, nkind
1669  iac = ikind + nkind*(kkind - 1)
1670  ibc = jkind + nkind*(kkind - 1)
1671  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1672  IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1673  CALL get_alist(sap_int(iac), alist_ac, iatom)
1674  CALL get_alist(sap_int(ibc), alist_bc, jatom)
1675  IF (.NOT. ASSOCIATED(alist_ac)) cycle
1676  IF (.NOT. ASSOCIATED(alist_bc)) cycle
1677  DO kac = 1, alist_ac%nclist
1678  DO kbc = 1, alist_bc%nclist
1679  IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1680  IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1681  IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1682 
1683  acint => alist_ac%clist(kac)%acint
1684  bcint => alist_bc%clist(kbc)%acint
1685  achint => alist_ac%clist(kac)%achint
1686  bchint => alist_bc%clist(kbc)%achint
1687  na = SIZE(acint, 1)
1688  np = SIZE(acint, 2)
1689  nb = SIZE(bcint, 1)
1690  ! Position of the pseudized atom
1691  r_ps = particle_set(alist_ac%clist(kac)%catom)%r
1692  r_b = refpoint
1693 
1694 !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1695 !$ CALL omp_set_lock(locks(hash))
1696  ! assemble integrals
1697  IF (iatom <= jatom) THEN
1698  blocks_mag(1)%block(1:na, 1:nb) = blocks_mag(1)%block(1:na, 1:nb) + &
1699  (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1700  matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) & ! R_y [V_nl, z]
1701  - (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1702  matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) ! - R_z [V_nl, y]
1703  blocks_mag(2)%block(1:na, 1:nb) = blocks_mag(2)%block(1:na, 1:nb) + &
1704  (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1705  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) & ! R_z [V_nl, x]
1706  - (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1707  matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) ! - R_x [V_nl, z]
1708  blocks_mag(3)%block(1:na, 1:nb) = blocks_mag(3)%block(1:na, 1:nb) + &
1709  (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1710  matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) & ! R_x [V_nl, y]
1711  - (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1712  matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) ! - R_y [V_nl, x]
1713  ELSE
1714  blocks_mag(1)%block(1:nb, 1:na) = blocks_mag(1)%block(1:nb, 1:na) + &
1715  (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1716  matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) & ! R_y [V_nl, z]
1717  - (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1718  matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) ! - R_z [V_nl, y]
1719  blocks_mag(2)%block(1:nb, 1:na) = blocks_mag(2)%block(1:nb, 1:na) + &
1720  (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1721  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) & ! R_z [V_nl, x]
1722  - (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1723  matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) ! - R_x [V_nl, z]
1724  blocks_mag(3)%block(1:nb, 1:na) = blocks_mag(3)%block(1:nb, 1:na) + &
1725  (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1726  matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) & ! R_x [V_nl, y]
1727  - (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1728  matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) ! - R_y [V_nl, x]
1729  END IF
1730 !$ CALL omp_unset_lock(locks(hash))
1731  EXIT ! We have found a match and there can be only one single match
1732  END IF
1733  END DO
1734  END DO
1735  END DO
1736  END IF
1737 
1738  DO ind = 1, 3
1739  NULLIFY (blocks_mag(ind)%block)
1740  END DO
1741  DEALLOCATE (blocks_mag)
1742  END DO
1743 
1744 !$OMP DO
1745 !$ DO lock_num = 1, nlock
1746 !$ call omp_destroy_lock(locks(lock_num))
1747 !$ END DO
1748 !$OMP END DO
1749 
1750 !$OMP SINGLE
1751 !$ DEALLOCATE (locks)
1752 !$OMP END SINGLE NOWAIT
1753 
1754 !$OMP END PARALLEL
1755 
1756  DEALLOCATE (basis_set)
1757  CALL release_sap_int(sap_int)
1758 
1759  CALL timestop(handle)
1760 
1761  END SUBROUTINE build_com_nl_mag
1762 
1763 ! **************************************************************************************************
1764 !> \brief Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs
1765 !> \param qs_kind_set ...
1766 !> \param sab_all ...
1767 !> \param sap_ppnl ...
1768 !> \param eps_ppnl ...
1769 !> \param particle_set ...
1770 !> \param matrix_rv ...
1771 !> \param ref_point ...
1772 !> \param cell ...
1773 !> \param direction_Or If set to true: calculate Vnl * r_delta
1774 !> Otherwise calculate r_delta * Vnl
1775 ! **************************************************************************************************
1776  SUBROUTINE build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, &
1777  matrix_rv, ref_point, cell, direction_Or)
1778 
1779  TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1780  POINTER :: qs_kind_set
1781  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1782  INTENT(IN), POINTER :: sab_all, sap_ppnl
1783  REAL(kind=dp), INTENT(IN) :: eps_ppnl
1784  TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1785  POINTER :: particle_set
1786  TYPE(dbcsr_p_type), DIMENSION(:, :), &
1787  INTENT(INOUT), OPTIONAL, POINTER :: matrix_rv
1788  REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
1789  TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1790  LOGICAL :: direction_or
1791 
1792  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_vnl_giao'
1793  INTEGER, PARAMETER :: i_1 = 1
1794 
1795  INTEGER :: delta, gamma, handle, i, iab, iac, &
1796  iatom, ibc, icol, ikind, irow, j, &
1797  jatom, jkind, kac, kbc, kkind, na, &
1798  natom, nb, nkind, np, order, slot
1799  INTEGER, DIMENSION(3) :: cell_b
1800  LOGICAL :: found, my_ref, ppnl_present
1801  REAL(kind=dp), DIMENSION(3) :: rab, rf
1802  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1803  TYPE(alist_type), POINTER :: alist_ac, alist_bc
1804  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1805  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1806  DIMENSION(:) :: basis_set
1807  TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1808  TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1809 
1810 !$ INTEGER(kind=omp_lock_kind), &
1811 !$ ALLOCATABLE, DIMENSION(:) :: locks
1812 !$ INTEGER :: lock_num, hash
1813 !$ INTEGER, PARAMETER :: nlock = 501
1814 
1815  ppnl_present = ASSOCIATED(sap_ppnl)
1816  IF (.NOT. ppnl_present) RETURN
1817 
1818  CALL timeset(routinen, handle)
1819 
1820  natom = SIZE(particle_set)
1821 
1822  my_ref = .false.
1823  IF (PRESENT(ref_point)) THEN
1824  cpassert(PRESENT(cell)) ! need cell as well if refpoint is provided
1825  rf = ref_point
1826  my_ref = .true.
1827  END IF
1828 
1829  nkind = SIZE(qs_kind_set)
1830 
1831  ! sap_int needs to be shared as multiple threads need to access this
1832  NULLIFY (sap_int)
1833  ALLOCATE (sap_int(nkind*nkind))
1834  DO i = 1, nkind*nkind
1835  NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1836  sap_int(i)%nalist = 0
1837  END DO
1838 
1839  order = 1
1840  IF (my_ref) THEN
1841  ! calculate integrals <a|x^n|p>
1842  CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=rf, &
1843  particle_set=particle_set, cell=cell)
1844  ELSE
1845  CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1846  END IF
1847 
1848  ! *** Set up a sorting index
1849  CALL sap_sort(sap_int)
1850 
1851  ALLOCATE (basis_set(nkind))
1852  DO ikind = 1, nkind
1853  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1854  IF (ASSOCIATED(orb_basis_set)) THEN
1855  basis_set(ikind)%gto_basis_set => orb_basis_set
1856  ELSE
1857  NULLIFY (basis_set(ikind)%gto_basis_set)
1858  END IF
1859  END DO
1860 
1861  CALL get_neighbor_list_set_p(neighbor_list_sets=sab_all)
1862  ! *** All integrals needed have been calculated and stored in sap_int
1863  ! *** We now calculate the commutator matrix elements
1864 
1865 !$OMP PARALLEL &
1866 !$OMP DEFAULT (NONE) &
1867 !$OMP SHARED (basis_set, matrix_rv, &
1868 !$OMP sap_int, nkind, eps_ppnl, locks, sab_all, &
1869 !$OMP particle_set, direction_Or) &
1870 !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1871 !$OMP iab, irow, icol, blocks_rv, &
1872 !$OMP found, iac, ibc, alist_ac, alist_bc, &
1873 !$OMP na, np, nb, kkind, kac, kbc, i, &
1874 !$OMP hash, natom, delta, gamma, achint, bchint, acint, bcint)
1875 
1876 !$OMP SINGLE
1877 !$ ALLOCATE (locks(nlock))
1878 !$OMP END SINGLE
1879 
1880 !$OMP DO
1881 !$ DO lock_num = 1, nlock
1882 !$ call omp_init_lock(locks(lock_num))
1883 !$ END DO
1884 !$OMP END DO
1885 
1886 !$OMP DO SCHEDULE(GUIDED)
1887 
1888  DO slot = 1, sab_all(1)%nl_size
1889 
1890  ikind = sab_all(1)%nlist_task(slot)%ikind
1891  jkind = sab_all(1)%nlist_task(slot)%jkind
1892  iatom = sab_all(1)%nlist_task(slot)%iatom
1893  jatom = sab_all(1)%nlist_task(slot)%jatom
1894  cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1895  rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1896 
1897  IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1898  IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1899  iab = ikind + nkind*(jkind - 1)
1900 
1901  irow = iatom
1902  icol = jatom
1903 
1904  ! allocate blocks
1905  ALLOCATE (blocks_rv(3, 3))
1906 
1907  ! get blocks
1908  DO i = 1, 3
1909  DO j = 1, 3
1910  CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, &
1911  blocks_rv(i, j)%block, found)
1912  blocks_rv(i, j)%block(:, :) = 0.0_dp
1913  cpassert(found)
1914  END DO
1915  END DO
1916 
1917  ! loop over all kinds for projector atom
1918  ! < iatom | katom > h < katom | jatom >
1919  DO kkind = 1, nkind
1920  iac = ikind + nkind*(kkind - 1)
1921  ibc = jkind + nkind*(kkind - 1)
1922  IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1923  IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1924  CALL get_alist(sap_int(iac), alist_ac, iatom)
1925  CALL get_alist(sap_int(ibc), alist_bc, jatom)
1926  IF (.NOT. ASSOCIATED(alist_ac)) cycle
1927  IF (.NOT. ASSOCIATED(alist_bc)) cycle
1928  DO kac = 1, alist_ac%nclist
1929  DO kbc = 1, alist_bc%nclist
1930  IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1931 
1932  IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1933  IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1934  acint => alist_ac%clist(kac)%acint
1935  bcint => alist_bc%clist(kbc)%acint
1936  achint => alist_ac%clist(kac)%achint
1937  bchint => alist_bc%clist(kbc)%achint
1938  na = SIZE(acint, 1)
1939  np = SIZE(acint, 2)
1940  nb = SIZE(bcint, 1)
1941 !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1942 !$ CALL omp_set_lock(locks(hash))
1943 
1944  !! The atom index is alist_ac%clist(kac)%catom
1945  ! The coordinate is particle_set(alist_ac%clist(kac)%catom)%r(:)
1946  IF (direction_or) THEN ! V * r_delta * (R^eta_gamma - R^nu_gamma)
1947  DO delta = 1, 3
1948  DO gamma = 1, 3
1949  blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1950  = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1951  matmul(achint(1:na, 1:np, i_1), transpose(bcint(1:nb, 1:np, delta + 1))) &
1952  *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1953  END DO
1954  END DO
1955  ELSE ! r_delta * V * (R^eta_gamma - R^nu_gamma)
1956  DO delta = 1, 3
1957  DO gamma = 1, 3
1958  blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1959  = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1960  matmul(achint(1:na, 1:np, delta + 1), transpose(bcint(1:nb, 1:np, i_1))) &
1961  *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1962  END DO
1963  END DO
1964  END IF
1965 
1966 !$ CALL omp_unset_lock(locks(hash))
1967  EXIT ! We have found a match and there can be only one single match
1968  END IF
1969  END DO
1970  END DO
1971  END DO
1972  DO delta = 1, 3
1973  DO gamma = 1, 3
1974  NULLIFY (blocks_rv(gamma, delta)%block)
1975  END DO
1976  END DO
1977  DEALLOCATE (blocks_rv)
1978  END DO
1979 
1980 !$OMP DO
1981 !$ DO lock_num = 1, nlock
1982 !$ call omp_destroy_lock(locks(lock_num))
1983 !$ END DO
1984 !$OMP END DO
1985 
1986 !$OMP SINGLE
1987 !$ DEALLOCATE (locks)
1988 !$OMP END SINGLE NOWAIT
1989 
1990 !$OMP END PARALLEL
1991 
1992  CALL release_sap_int(sap_int)
1993 
1994  DEALLOCATE (basis_set)
1995 
1996  CALL timestop(handle)
1997 
1998  END SUBROUTINE build_com_vnl_giao
1999 
2000 END MODULE commutator_rpnl
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition: ai_moments.F:17
ai_moments::moment
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab)
...
Definition: ai_moments.F:986
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
ai_overlap::overlap
subroutine, public overlap(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, rab, dab, sab, da_max_set, return_derivatives, s, lds, sdab, pab, force_a)
Purpose: Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions...
Definition: ai_overlap.F:73
basis_set_types
Definition: basis_set_types.F:15
block_p_types
collect pointers to a block of reals
Definition: block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: commutator_rpnl.F:14
commutator_rpnl::build_com_vnl_giao
subroutine, public build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_rv, ref_point, cell, direction_Or)
Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs.
Definition: commutator_rpnl.F:1778
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition: commutator_rpnl.F:439
commutator_rpnl::build_com_nl_mag
subroutine, public build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
Definition: commutator_rpnl.F:1538
commutator_rpnl::build_com_rpnl
subroutine, public build_com_rpnl(matrix_rv, qs_kind_set, sab_orb, sap_ppnl, eps_ppnl)
...
Definition: commutator_rpnl.F:76
external_potential_types
Definition of the atomic potential types.
Definition: external_potential_types.F:14
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition: orbital_pointers.F:253
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition: orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_neighbor_list_set_p
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Definition: qs_neighbor_list_types.F:993
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
sap_kind_types
General overlap type integrals containers.
Definition: sap_kind_types.F:12
sap_kind_types::build_sap_ints
subroutine, public build_sap_ints(sap_int, sap_ppnl, qs_kind_set, nder, moment_mode, refpoint, particle_set, cell, pseudoatom)
Calculate overlap and optionally momenta <a|x^n|p> between GTOs and nl. pseudo potential projectors a...
Definition: sap_kind_types.F:282
sap_kind_types::release_sap_int
subroutine, public release_sap_int(sap_int)
...
Definition: sap_kind_types.F:86
sap_kind_types::sap_sort
subroutine, public sap_sort(sap_int)
...
Definition: sap_kind_types.F:247
sap_kind_types::get_alist
subroutine, public get_alist(sap_int, alist, atom)
...
Definition: sap_kind_types.F:140