(git:ed6f26b)
Loading...
Searching...
No Matches
commutator_rpnl.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian
9!> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
10!> \par History
11!> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01]
12!> - full rewrite [jhu, 2009-01-23]
13! **************************************************************************************************
14MODULE commutator_rpnl
15 USE 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 cp_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
65CONTAINS
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
2000END 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:206
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_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_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_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
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
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, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, 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_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
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, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
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
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
block_p_types::block_p_type
Definition block_p_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
external_potential_types::gth_potential_p_type
Definition external_potential_types.F:227
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
external_potential_types::sgp_potential_p_type
Definition external_potential_types.F:235
external_potential_types::sgp_potential_type
Definition external_potential_types.F:172
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
sap_kind_types::alist_type
Definition sap_kind_types.F:58
sap_kind_types::clist_type
Definition sap_kind_types.F:43
sap_kind_types::sap_int_type
Definition sap_kind_types.F:64