(git:8917686)
Loading...
Searching...
No Matches
commutator_rpnl.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8!> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian
9!> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
10!> \par History
11!> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01]
12!> - full rewrite [jhu, 2009-01-23]
13! **************************************************************************************************
14MODULE commutator_rpnl
15 USE 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 ELSE IF (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 > 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, lock_num, &
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 IF (iatom <= jatom) THEN
872 ! yzV
873 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
874 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
875 ! -yVz
876 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
877 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4)))
878 ! -zyV
879 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
880 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
881 ! zVy
882 blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
883 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 3)))
884 ELSE
885 ! yzV
886 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
887 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
888 ! -yVz
889 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
890 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4)))
891 ! -zyV
892 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
893 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
894 ! zVy
895 blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
896 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 3)))
897 END IF
898
899 ! y-component (z [x,Vnl] - x [z, Vnl])
900 IF (iatom <= jatom) THEN
901 ! zxV
902 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
903 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
904 ! -zVx
905 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
906 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 2)))
907 ! -xzV
908 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
909 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
910 ! xVz
911 blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
912 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4)))
913 ELSE
914 ! zxV
915 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
916 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
917 ! -zVx
918 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
919 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 2)))
920 ! -xzV
921 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
922 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
923 ! xVz
924 blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
925 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4)))
926 END IF
927
928 ! z-component (x [y,Vnl] - y [x, Vnl])
929 IF (iatom <= jatom) THEN
930 ! xyV
931 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
932 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
933 ! -xVy
934 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
935 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3)))
936 ! -yxV
937 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
938 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
939 ! zVx
940 blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
941 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 2)))
942 ELSE
943 ! xyV
944 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
945 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
946 ! -xVy
947 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
948 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3)))
949 ! -yxV
950 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
951 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
952 ! zVx
953 blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
954 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 2)))
955 END IF
956 END IF
957
958 IF (my_rrv) THEN
959 ! r_alpha * r_beta * Vnl
960 IF (iatom <= jatom) THEN
961 ! xxV
962 blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) + &
963 matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1)))
964 ! xyV
965 blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) + &
966 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
967 ! xzV
968 blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) + &
969 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
970 ! yyV
971 blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) + &
972 matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1)))
973 ! yzV
974 blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) + &
975 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
976 ! zzV
977 blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) + &
978 matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1)))
979 ELSE
980 ! xxV
981 blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) + &
982 matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1)))
983 ! xyV
984 blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) + &
985 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
986 ! xzV
987 blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) + &
988 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
989 ! yyV
990 blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) + &
991 matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1)))
992 ! yzV
993 blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) + &
994 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
995 ! zzV
996 blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) + &
997 matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1)))
998 END IF
999
1000 ! - Vnl * r_alpha * r_beta
1001 IF (iatom <= jatom) THEN
1002 ! -Vxx
1003 blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) - &
1004 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5)))
1005 ! -Vxy
1006 blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) - &
1007 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6)))
1008 ! -Vxz
1009 blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) - &
1010 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7)))
1011 ! -Vyy
1012 blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) - &
1013 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8)))
1014 ! -Vyz
1015 blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) - &
1016 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9)))
1017 ! -Vzz
1018 blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) - &
1019 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10)))
1020 ELSE
1021 ! -Vxx
1022 blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) - &
1023 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5)))
1024 ! -Vxy
1025 blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) - &
1026 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6)))
1027 ! -Vxz
1028 blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) - &
1029 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7)))
1030 ! -Vyy
1031 blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) - &
1032 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8)))
1033 ! -Vyz
1034 blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) - &
1035 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9)))
1036 ! -Vzz
1037 blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) - &
1038 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10)))
1039 END IF
1040 END IF
1041
1042 IF (my_rvr) THEN
1043 ! r_alpha * Vnl * r_beta
1044 IF (iatom <= jatom) THEN
1045 ! xVx
1046 blocks_rvr(1)%block(1:na, 1:nb) = blocks_rvr(1)%block(1:na, 1:nb) + &
1047 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 2)))
1048 ! xVy
1049 blocks_rvr(2)%block(1:na, 1:nb) = blocks_rvr(2)%block(1:na, 1:nb) + &
1050 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 3)))
1051 ! xVz
1052 blocks_rvr(3)%block(1:na, 1:nb) = blocks_rvr(3)%block(1:na, 1:nb) + &
1053 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 4)))
1054 ! yVy
1055 blocks_rvr(4)%block(1:na, 1:nb) = blocks_rvr(4)%block(1:na, 1:nb) + &
1056 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 3)))
1057 ! yVz
1058 blocks_rvr(5)%block(1:na, 1:nb) = blocks_rvr(5)%block(1:na, 1:nb) + &
1059 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 4)))
1060 ! zVz
1061 blocks_rvr(6)%block(1:na, 1:nb) = blocks_rvr(6)%block(1:na, 1:nb) + &
1062 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 4)))
1063 ELSE
1064 ! xVx
1065 blocks_rvr(1)%block(1:nb, 1:na) = blocks_rvr(1)%block(1:nb, 1:na) + &
1066 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 2)))
1067 ! xVy
1068 blocks_rvr(2)%block(1:nb, 1:na) = blocks_rvr(2)%block(1:nb, 1:na) + &
1069 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 3)))
1070 ! xVz
1071 blocks_rvr(3)%block(1:nb, 1:na) = blocks_rvr(3)%block(1:nb, 1:na) + &
1072 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 4)))
1073 ! yVy
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)))
1076 ! yVz
1077 blocks_rvr(5)%block(1:nb, 1:na) = blocks_rvr(5)%block(1:nb, 1:na) + &
1078 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 4)))
1079 ! zVz
1080 blocks_rvr(6)%block(1:nb, 1:na) = blocks_rvr(6)%block(1:nb, 1:na) + &
1081 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 4)))
1082 END IF
1083 END IF
1084
1085 IF (my_rrv_vrr) THEN
1086 ! r_alpha * r_beta * Vnl
1087 IF (iatom <= jatom) THEN
1088 ! xxV
1089 blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
1090 matmul(achint(1:na, 1:np, 5), transpose(bcint(1:nb, 1:np, 1)))
1091 ! xyV
1092 blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
1093 matmul(achint(1:na, 1:np, 6), transpose(bcint(1:nb, 1:np, 1)))
1094 ! xzV
1095 blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
1096 matmul(achint(1:na, 1:np, 7), transpose(bcint(1:nb, 1:np, 1)))
1097 ! yyV
1098 blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
1099 matmul(achint(1:na, 1:np, 8), transpose(bcint(1:nb, 1:np, 1)))
1100 ! yzV
1101 blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
1102 matmul(achint(1:na, 1:np, 9), transpose(bcint(1:nb, 1:np, 1)))
1103 ! zzV
1104 blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
1105 matmul(achint(1:na, 1:np, 10), transpose(bcint(1:nb, 1:np, 1)))
1106 ELSE
1107 ! xxV
1108 blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
1109 matmul(bchint(1:nb, 1:np, 5), transpose(acint(1:na, 1:np, 1)))
1110 ! xyV
1111 blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
1112 matmul(bchint(1:nb, 1:np, 6), transpose(acint(1:na, 1:np, 1)))
1113 ! xzV
1114 blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
1115 matmul(bchint(1:nb, 1:np, 7), transpose(acint(1:na, 1:np, 1)))
1116 ! yyV
1117 blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
1118 matmul(bchint(1:nb, 1:np, 8), transpose(acint(1:na, 1:np, 1)))
1119 ! yzV
1120 blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
1121 matmul(bchint(1:nb, 1:np, 9), transpose(acint(1:na, 1:np, 1)))
1122 ! zzV
1123 blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
1124 matmul(bchint(1:nb, 1:np, 10), transpose(acint(1:na, 1:np, 1)))
1125 END IF
1126 ! + Vnl * r_alpha * r_beta
1127 IF (iatom <= jatom) THEN
1128 ! +Vxx
1129 blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
1130 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 5)))
1131 ! +Vxy
1132 blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
1133 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 6)))
1134 ! +Vxz
1135 blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
1136 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 7)))
1137 ! +Vyy
1138 blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
1139 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 8)))
1140 ! +Vyz
1141 blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
1142 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 9)))
1143 ! +Vzz
1144 blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
1145 matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 10)))
1146 ELSE
1147 ! +Vxx
1148 blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
1149 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 5)))
1150 ! +Vxy
1151 blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
1152 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 6)))
1153 ! +Vxz
1154 blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
1155 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 7)))
1156 ! +Vyy
1157 blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
1158 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 8)))
1159 ! +Vyz
1160 blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
1161 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 9)))
1162 ! +Vzz
1163 blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
1164 matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 10)))
1165 END IF
1166 END IF
1167
1168 ! The indices are stored in i_1, i_x, ..., i_zzz
1169
1170 ! TODO: is this set to zero before?
1171 IF (my_r_rxvr) THEN
1172 ! beta = 1
1173 ! matrix_r_rxvr(x, x) = x * y * V_nl * z - x * z * V_nl * y
1174 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
1175 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) + &
1176 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
1177 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
1178 blocks_r_rxvr(1, 1)%block(1:na, 1:nb) - &
1179 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
1180
1181 ! matrix_r_rxvr(y, x) = x * z * V_nl * x - x * x * V_nl * z
1182 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
1183 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) + &
1184 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_x)))
1185 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
1186 blocks_r_rxvr(2, 1)%block(1:na, 1:nb) - &
1187 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
1188
1189 ! matrix_r_rxvr(z, x) = x * x * V_nl * y - x * y * V_nl * x
1190 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
1191 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) + &
1192 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
1193 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
1194 blocks_r_rxvr(3, 1)%block(1:na, 1:nb) - &
1195 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_x)))
1196
1197 ! beta = 2
1198 ! matrix_r_rxvr(x, y) = y * y * V_nl * z - y * z * V_nl * y
1199 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
1200 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) + &
1201 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
1202 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
1203 blocks_r_rxvr(1, 2)%block(1:na, 1:nb) - &
1204 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_y)))
1205
1206 ! matrix_r_rxvr(y, y) = y * z * V_nl * x - y * x * V_nl * z
1207 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
1208 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) + &
1209 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
1210 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
1211 blocks_r_rxvr(2, 2)%block(1:na, 1:nb) - &
1212 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
1213
1214 ! matrix_r_rxvr(z, y) = y * x * V_nl * y - y * y * V_nl * x
1215 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
1216 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) + &
1217 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_y)))
1218 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
1219 blocks_r_rxvr(3, 2)%block(1:na, 1:nb) - &
1220 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
1221
1222 ! beta = 3
1223 ! matrix_r_rxvr(x, z) = z * y * V_nl * z - z * z * V_nl * y
1224 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
1225 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) + &
1226 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_z)))
1227 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
1228 blocks_r_rxvr(1, 3)%block(1:na, 1:nb) - &
1229 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
1230
1231 ! matrix_r_rxvr(y, z) = z * z * V_nl * x - z * x * V_nl * z
1232 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
1233 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) + &
1234 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
1235 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
1236 blocks_r_rxvr(2, 3)%block(1:na, 1:nb) - &
1237 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_z)))
1238
1239 ! matrix_r_rxvr(z, z) = z * x * V_nl * y - z * y * V_nl * x
1240 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
1241 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) + &
1242 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
1243 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
1244 blocks_r_rxvr(3, 3)%block(1:na, 1:nb) - &
1245 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
1246
1247 END IF ! my_r_rxvr
1248
1249 ! The indices are stored in i_1, i_x, ..., i_zzz
1250 ! This will put into blocks_rxvr_r
1251 ! matrix_rxvr_r(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
1252 ! r_gamma * V_nl * r_delta * r_beta
1253 IF (my_rxvr_r) THEN
1254 ! beta = 1
1255 ! matrix_rxvr_r(x, x) = yV zx - zV yx
1256 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
1257 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) + &
1258 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
1259 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
1260 blocks_rxvr_r(1, 1)%block(1:na, 1:nb) - &
1261 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
1262
1263 ! matrix_rxvr_r(y, x) = zV xx - xV zx
1264 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
1265 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) + &
1266 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
1267 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
1268 blocks_rxvr_r(2, 1)%block(1:na, 1:nb) - &
1269 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zx)))
1270
1271 ! matrix_rxvr_r(z, x) = xV yx - yV xx
1272 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
1273 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) + &
1274 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yx)))
1275 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
1276 blocks_rxvr_r(3, 1)%block(1:na, 1:nb) - &
1277 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
1278
1279 ! beta = 2
1280 ! matrix_rxvr_r(x, y) = yV zy - zV yy
1281 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
1282 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) + &
1283 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zy)))
1284 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
1285 blocks_rxvr_r(1, 2)%block(1:na, 1:nb) - &
1286 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
1287
1288 ! matrix_rxvr_r(y, y) = zV xy - xV zy
1289 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
1290 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) + &
1291 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
1292 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
1293 blocks_rxvr_r(2, 2)%block(1:na, 1:nb) - &
1294 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
1295
1296 ! matrix_rxvr_r(z, y) = xV yy - yV xy
1297 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
1298 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) + &
1299 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
1300 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
1301 blocks_rxvr_r(3, 2)%block(1:na, 1:nb) - &
1302 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xy)))
1303
1304 ! beta = 3
1305 ! matrix_rxvr_r(x, z) = yV zz - zV yz
1306 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
1307 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) + &
1308 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
1309 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
1310 blocks_rxvr_r(1, 3)%block(1:na, 1:nb) - &
1311 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yz)))
1312
1313 ! matrix_rxvr_r(y, z) = zV xz - xV zz
1314 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
1315 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) + &
1316 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xz)))
1317 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
1318 blocks_rxvr_r(2, 3)%block(1:na, 1:nb) - &
1319 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
1320
1321 ! matrix_rxvr_r(z, z) = xV yz - yV xz
1322 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
1323 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) + &
1324 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
1325 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
1326 blocks_rxvr_r(3, 3)%block(1:na, 1:nb) - &
1327 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
1328
1329 END IF ! my_rxvr_r
1330
1331 ! matrix_r_doublecom(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
1332 ! gamma V^pseudoatom beta delta - gamma beta V^pseudoatom delta
1333
1334 IF (my_r_doublecom) THEN
1335 ! beta = 1
1336 ! matrix_r_doublecom(x, x) = yV xz - zV xy - yxV z + zxV y
1337 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1338 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
1339 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xz)))
1340 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1341 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
1342 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xy)))
1343 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1344 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
1345 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_z)))
1346 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1347 blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
1348 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_y)))
1349
1350 ! matrix_r_doublecom(y, x) = zV xx - xV xz - zxV x + xxV z
1351 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1352 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
1353 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_xx)))
1354 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1355 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
1356 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xz)))
1357 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1358 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
1359 matmul(achint(1:na, 1:np, i_zx), transpose(bcint(1:nb, 1:np, i_x)))
1360 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1361 blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
1362 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_z)))
1363
1364 ! matrix_r_doublecom(z, x) = xV xy - yV xx - xxV y + yxV x
1365 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1366 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1367 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_xy)))
1368 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1369 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1370 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_xx)))
1371 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1372 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1373 matmul(achint(1:na, 1:np, i_xx), transpose(bcint(1:nb, 1:np, i_y)))
1374 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1375 blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1376 matmul(achint(1:na, 1:np, i_yx), transpose(bcint(1:nb, 1:np, i_x)))
1377
1378 ! beta = 2
1379 ! matrix_r_doublecom(x, y) = yV yz - zV yy - yyV z + zyV y
1380 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1381 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1382 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yz)))
1383 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1384 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1385 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yy)))
1386 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1387 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1388 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_z)))
1389 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1390 blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1391 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_y)))
1392
1393 ! matrix_r_doublecom(y, y) = zV yx - xV yz - zyV x + xyV z
1394 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1395 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1396 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_yx)))
1397 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1398 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1399 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yz)))
1400 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1401 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1402 matmul(achint(1:na, 1:np, i_zy), transpose(bcint(1:nb, 1:np, i_x)))
1403 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1404 blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1405 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_z)))
1406
1407 ! matrix_r_doublecom(z, y) = xV yy - yV yx - xyV y + yyV x
1408 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1409 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1410 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_yy)))
1411 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1412 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1413 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_yx)))
1414 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1415 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1416 matmul(achint(1:na, 1:np, i_xy), transpose(bcint(1:nb, 1:np, i_y)))
1417 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1418 blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1419 matmul(achint(1:na, 1:np, i_yy), transpose(bcint(1:nb, 1:np, i_x)))
1420
1421 ! beta = 3
1422 ! matrix_r_doublecom(x, z) = yV zz - zV zy - yzV z + zzV y
1423 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1424 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1425 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zz)))
1426 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1427 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1428 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zy)))
1429 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1430 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1431 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_z)))
1432 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1433 blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1434 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_y)))
1435
1436 ! matrix_r_doublecom(y, z) = zV zx - xV zz - zzV x + xzV z
1437 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1438 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1439 matmul(achint(1:na, 1:np, i_z), transpose(bcint(1:nb, 1:np, i_zx)))
1440 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1441 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1442 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zz)))
1443 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1444 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1445 matmul(achint(1:na, 1:np, i_zz), transpose(bcint(1:nb, 1:np, i_x)))
1446 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1447 blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1448 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_z)))
1449
1450 ! matrix_r_doublecom(z, z) = xV zy - yV zx - xzV y + yzV x
1451 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1452 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1453 matmul(achint(1:na, 1:np, i_x), transpose(bcint(1:nb, 1:np, i_zy)))
1454 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1455 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1456 matmul(achint(1:na, 1:np, i_y), transpose(bcint(1:nb, 1:np, i_zx)))
1457 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1458 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1459 matmul(achint(1:na, 1:np, i_xz), transpose(bcint(1:nb, 1:np, i_y)))
1460 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1461 blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1462 matmul(achint(1:na, 1:np, i_yz), transpose(bcint(1:nb, 1:np, i_x)))
1463
1464 END IF ! my_r_doublecom
1465!$ CALL omp_unset_lock(locks(hash))
1466 EXIT ! We have found a match and there can be only one single match
1467 END IF
1468 END DO
1469 END DO
1470 END DO
1471 END IF
1472 IF (my_rv) THEN
1473 DO ind = 1, 3
1474 NULLIFY (blocks_rv(ind)%block)
1475 END DO
1476 DEALLOCATE (blocks_rv)
1477 END IF
1478 IF (my_rxrv) THEN
1479 DO ind = 1, 3
1480 NULLIFY (blocks_rxrv(ind)%block)
1481 END DO
1482 DEALLOCATE (blocks_rxrv)
1483 END IF
1484 IF (my_rrv) THEN
1485 DO ind = 1, 6
1486 NULLIFY (blocks_rrv(ind)%block)
1487 END DO
1488 DEALLOCATE (blocks_rrv)
1489 END IF
1490 IF (my_rvr) THEN
1491 DO ind = 1, 6
1492 NULLIFY (blocks_rvr(ind)%block)
1493 END DO
1494 DEALLOCATE (blocks_rvr)
1495 END IF
1496 IF (my_rrv_vrr) THEN
1497 DO ind = 1, 6
1498 NULLIFY (blocks_rrv_vrr(ind)%block)
1499 END DO
1500 DEALLOCATE (blocks_rrv_vrr)
1501 END IF
1502 IF (my_r_rxvr) THEN
1503 DO ind = 1, 3
1504 DO ind2 = 1, 3
1505 NULLIFY (blocks_r_rxvr(ind, ind2)%block)
1506 END DO
1507 END DO
1508 DEALLOCATE (blocks_r_rxvr)
1509 END IF
1510 IF (my_rxvr_r) THEN
1511 DO ind = 1, 3
1512 DO ind2 = 1, 3
1513 NULLIFY (blocks_rxvr_r(ind, ind2)%block)
1514 END DO
1515 END DO
1516 DEALLOCATE (blocks_rxvr_r)
1517 END IF
1518 IF (my_r_doublecom) THEN
1519 DO ind = 1, 3
1520 DO ind2 = 1, 3
1521 NULLIFY (blocks_r_doublecom(ind, ind2)%block)
1522 END DO
1523 END DO
1524 DEALLOCATE (blocks_r_doublecom)
1525 END IF
1526 END DO
1527
1528!$OMP DO
1529!$ DO lock_num = 1, nlock
1530!$ call omp_destroy_lock(locks(lock_num))
1531!$ END DO
1532!$OMP END DO
1533
1534!$OMP SINGLE
1535!$ DEALLOCATE (locks)
1536!$OMP END SINGLE NOWAIT
1537
1538!$OMP END PARALLEL
1539
1540 CALL release_sap_int(sap_int)
1541
1542 DEALLOCATE (basis_set)
1543
1544 CALL timestop(handle)
1545
1546 END SUBROUTINE build_com_mom_nl
1547
1548! **************************************************************************************************
1549!> \brief calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
1550!> \param qs_kind_set ...
1551!> \param sab_all ...
1552!> \param sap_ppnl ...
1553!> \param eps_ppnl ...
1554!> \param particle_set ...
1555!> \param matrix_mag_nl ...
1556!> \param refpoint ...
1557!> \param cell ...
1558! **************************************************************************************************
1559 SUBROUTINE build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
1560
1561 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1562 POINTER :: qs_kind_set
1563 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1564 INTENT(IN), POINTER :: sab_all, sap_ppnl
1565 REAL(kind=dp), INTENT(IN) :: eps_ppnl
1566 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1567 POINTER :: particle_set
1568 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1569 POINTER :: matrix_mag_nl
1570 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
1571 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1572
1573 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_nl_mag'
1574
1575 INTEGER :: handle, iab, iac, iatom, ibc, icol, &
1576 ikind, ind, irow, jatom, jkind, kac, &
1577 kbc, kkind, na, natom, nb, nkind, np, &
1578 order, slot
1579 INTEGER, DIMENSION(3) :: cell_b
1580 LOGICAL :: found, go, my_ref, ppnl_present
1581 REAL(kind=dp), DIMENSION(3) :: r_b, r_ps, rab
1582 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1583 TYPE(alist_type), POINTER :: alist_ac, alist_bc
1584 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1585 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1586 DIMENSION(:) :: basis_set
1587 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1588 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1589
1590!$ INTEGER(kind=omp_lock_kind), &
1591!$ ALLOCATABLE, DIMENSION(:) :: locks
1592!$ INTEGER :: lock_num, hash
1593!$ INTEGER, PARAMETER :: nlock = 501
1594
1595 ppnl_present = ASSOCIATED(sap_ppnl)
1596 IF (.NOT. ppnl_present) RETURN
1597
1598 CALL timeset(routinen, handle)
1599
1600 my_ref = .false.
1601 IF (PRESENT(refpoint)) THEN
1602 my_ref = .true.
1603 cpassert(PRESENT(cell))
1604 END IF
1605
1606 natom = SIZE(particle_set)
1607 nkind = SIZE(qs_kind_set)
1608
1609 ! allocate integral storage
1610 NULLIFY (sap_int)
1611 ALLOCATE (sap_int(nkind*nkind))
1612 DO ind = 1, nkind*nkind
1613 NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1614 sap_int(ind)%nalist = 0
1615 END DO
1616
1617 ! build integrals over GTO + projector functions, refpoint actually
1618 order = 1 ! only need first moments (x, y, z)
1619 ! refpoint actually does not matter in this case, i. e. (order = 1 .and. commutator)
1620 IF (my_ref) THEN
1621 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=refpoint, &
1622 particle_set=particle_set, cell=cell)
1623 ELSE
1624 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1625 END IF
1626
1627 CALL sap_sort(sap_int)
1628
1629 ! get access to basis sets
1630 ALLOCATE (basis_set(nkind))
1631 DO ikind = 1, nkind
1632 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1633 IF (ASSOCIATED(orb_basis_set)) THEN
1634 basis_set(ikind)%gto_basis_set => orb_basis_set
1635 ELSE
1636 NULLIFY (basis_set(ikind)%gto_basis_set)
1637 END IF
1638 END DO
1639
1640!$OMP PARALLEL &
1641!$OMP DEFAULT (NONE) &
1642!$OMP SHARED (basis_set, matrix_mag_nl, sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
1643!$OMP particle_set, my_ref, refpoint) &
1644!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, lock_num, &
1645!$OMP iab, irow, icol, blocks_mag, r_ps, r_b, go, hash, &
1646!$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
1647!$OMP achint, bchint, na, np, nb, kkind, kac, kbc)
1648
1649!$OMP SINGLE
1650!$ ALLOCATE (locks(nlock))
1651!$OMP END SINGLE
1652
1653!$OMP DO
1654!$ DO lock_num = 1, nlock
1655!$ call omp_init_lock(locks(lock_num))
1656!$ END DO
1657!$OMP END DO
1658
1659!$OMP DO SCHEDULE(GUIDED)
1660 DO slot = 1, sab_all(1)%nl_size
1661 ! get indices
1662 ikind = sab_all(1)%nlist_task(slot)%ikind
1663 jkind = sab_all(1)%nlist_task(slot)%jkind
1664 iatom = sab_all(1)%nlist_task(slot)%iatom
1665 jatom = sab_all(1)%nlist_task(slot)%jatom
1666 cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1667 rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1668
1669 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1670 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1671 iab = ikind + nkind*(jkind - 1)
1672
1673 IF (iatom <= jatom) THEN
1674 irow = iatom
1675 icol = jatom
1676 ELSE
1677 irow = jatom
1678 icol = iatom
1679 END IF
1680
1681 ! get blocks
1682 ALLOCATE (blocks_mag(3))
1683 DO ind = 1, 3
1684 CALL dbcsr_get_block_p(matrix_mag_nl(ind)%matrix, irow, icol, blocks_mag(ind)%block, found)
1685 END DO
1686
1687 go = (ASSOCIATED(blocks_mag(1)%block) .AND. ASSOCIATED(blocks_mag(2)%block) .AND. ASSOCIATED(blocks_mag(3)%block))
1688
1689 IF (go) THEN
1690 DO kkind = 1, nkind
1691 iac = ikind + nkind*(kkind - 1)
1692 ibc = jkind + nkind*(kkind - 1)
1693 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1694 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1695 CALL get_alist(sap_int(iac), alist_ac, iatom)
1696 CALL get_alist(sap_int(ibc), alist_bc, jatom)
1697 IF (.NOT. ASSOCIATED(alist_ac)) cycle
1698 IF (.NOT. ASSOCIATED(alist_bc)) cycle
1699 DO kac = 1, alist_ac%nclist
1700 DO kbc = 1, alist_bc%nclist
1701 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1702 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1703 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1704
1705 acint => alist_ac%clist(kac)%acint
1706 bcint => alist_bc%clist(kbc)%acint
1707 achint => alist_ac%clist(kac)%achint
1708 bchint => alist_bc%clist(kbc)%achint
1709 na = SIZE(acint, 1)
1710 np = SIZE(acint, 2)
1711 nb = SIZE(bcint, 1)
1712 ! Position of the pseudized atom
1713 r_ps = particle_set(alist_ac%clist(kac)%catom)%r
1714 r_b = refpoint
1715
1716!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1717!$ CALL omp_set_lock(locks(hash))
1718 ! assemble integrals
1719 IF (iatom <= jatom) THEN
1720 blocks_mag(1)%block(1:na, 1:nb) = blocks_mag(1)%block(1:na, 1:nb) + &
1721 (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1722 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) & ! R_y [V_nl, z]
1723 - (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1724 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) ! - R_z [V_nl, y]
1725 blocks_mag(2)%block(1:na, 1:nb) = blocks_mag(2)%block(1:na, 1:nb) + &
1726 (r_ps(3) - r_b(3))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1727 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) & ! R_z [V_nl, x]
1728 - (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 4))) - &
1729 matmul(achint(1:na, 1:np, 4), transpose(bcint(1:nb, 1:np, 1)))) ! - R_x [V_nl, z]
1730 blocks_mag(3)%block(1:na, 1:nb) = blocks_mag(3)%block(1:na, 1:nb) + &
1731 (r_ps(1) - r_b(1))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 3))) - &
1732 matmul(achint(1:na, 1:np, 3), transpose(bcint(1:nb, 1:np, 1)))) & ! R_x [V_nl, y]
1733 - (r_ps(2) - r_b(2))*(matmul(achint(1:na, 1:np, 1), transpose(bcint(1:nb, 1:np, 2))) - &
1734 matmul(achint(1:na, 1:np, 2), transpose(bcint(1:nb, 1:np, 1)))) ! - R_y [V_nl, x]
1735 ELSE
1736 blocks_mag(1)%block(1:nb, 1:na) = blocks_mag(1)%block(1:nb, 1:na) + &
1737 (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1738 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) & ! R_y [V_nl, z]
1739 - (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1740 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) ! - R_z [V_nl, y]
1741 blocks_mag(2)%block(1:nb, 1:na) = blocks_mag(2)%block(1:nb, 1:na) + &
1742 (r_ps(3) - r_b(3))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1743 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) & ! R_z [V_nl, x]
1744 - (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 4))) - &
1745 matmul(bchint(1:nb, 1:np, 4), transpose(acint(1:na, 1:np, 1)))) ! - R_x [V_nl, z]
1746 blocks_mag(3)%block(1:nb, 1:na) = blocks_mag(3)%block(1:nb, 1:na) + &
1747 (r_ps(1) - r_b(1))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 3))) - &
1748 matmul(bchint(1:nb, 1:np, 3), transpose(acint(1:na, 1:np, 1)))) & ! R_x [V_nl, y]
1749 - (r_ps(2) - r_b(2))*(matmul(bchint(1:nb, 1:np, 1), transpose(acint(1:na, 1:np, 2))) - &
1750 matmul(bchint(1:nb, 1:np, 2), transpose(acint(1:na, 1:np, 1)))) ! - R_y [V_nl, x]
1751 END IF
1752!$ CALL omp_unset_lock(locks(hash))
1753 EXIT ! We have found a match and there can be only one single match
1754 END IF
1755 END DO
1756 END DO
1757 END DO
1758 END IF
1759
1760 DO ind = 1, 3
1761 NULLIFY (blocks_mag(ind)%block)
1762 END DO
1763 DEALLOCATE (blocks_mag)
1764 END DO
1765
1766!$OMP DO
1767!$ DO lock_num = 1, nlock
1768!$ call omp_destroy_lock(locks(lock_num))
1769!$ END DO
1770!$OMP END DO
1771
1772!$OMP SINGLE
1773!$ DEALLOCATE (locks)
1774!$OMP END SINGLE NOWAIT
1775
1776!$OMP END PARALLEL
1777
1778 DEALLOCATE (basis_set)
1779 CALL release_sap_int(sap_int)
1780
1781 CALL timestop(handle)
1782
1783 END SUBROUTINE build_com_nl_mag
1784
1785! **************************************************************************************************
1786!> \brief Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs
1787!> \param qs_kind_set ...
1788!> \param sab_all ...
1789!> \param sap_ppnl ...
1790!> \param eps_ppnl ...
1791!> \param particle_set ...
1792!> \param matrix_rv ...
1793!> \param ref_point ...
1794!> \param cell ...
1795!> \param direction_Or If set to true: calculate Vnl * r_delta
1796!> Otherwise calculate r_delta * Vnl
1797! **************************************************************************************************
1798 SUBROUTINE build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, &
1799 matrix_rv, ref_point, cell, direction_Or)
1800
1801 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1802 POINTER :: qs_kind_set
1803 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1804 INTENT(IN), POINTER :: sab_all, sap_ppnl
1805 REAL(kind=dp), INTENT(IN) :: eps_ppnl
1806 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1807 POINTER :: particle_set
1808 TYPE(dbcsr_p_type), DIMENSION(:, :), &
1809 INTENT(INOUT), OPTIONAL, POINTER :: matrix_rv
1810 REAL(kind=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
1811 TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1812 LOGICAL :: direction_or
1813
1814 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_com_vnl_giao'
1815 INTEGER, PARAMETER :: i_1 = 1
1816
1817 INTEGER :: delta, gamma, handle, i, iab, iac, &
1818 iatom, ibc, icol, ikind, irow, j, &
1819 jatom, jkind, kac, kbc, kkind, na, &
1820 natom, nb, nkind, np, order, slot
1821 INTEGER, DIMENSION(3) :: cell_b
1822 LOGICAL :: found, my_ref, ppnl_present
1823 REAL(kind=dp), DIMENSION(3) :: rab, rf
1824 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1825 TYPE(alist_type), POINTER :: alist_ac, alist_bc
1826 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1827 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1828 DIMENSION(:) :: basis_set
1829 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1830 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1831
1832!$ INTEGER(kind=omp_lock_kind), &
1833!$ ALLOCATABLE, DIMENSION(:) :: locks
1834!$ INTEGER :: lock_num, hash
1835!$ INTEGER, PARAMETER :: nlock = 501
1836
1837 ppnl_present = ASSOCIATED(sap_ppnl)
1838 IF (.NOT. ppnl_present) RETURN
1839
1840 CALL timeset(routinen, handle)
1841
1842 natom = SIZE(particle_set)
1843
1844 my_ref = .false.
1845 IF (PRESENT(ref_point)) THEN
1846 cpassert(PRESENT(cell)) ! need cell as well if refpoint is provided
1847 rf = ref_point
1848 my_ref = .true.
1849 END IF
1850
1851 nkind = SIZE(qs_kind_set)
1852
1853 ! sap_int needs to be shared as multiple threads need to access this
1854 NULLIFY (sap_int)
1855 ALLOCATE (sap_int(nkind*nkind))
1856 DO i = 1, nkind*nkind
1857 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1858 sap_int(i)%nalist = 0
1859 END DO
1860
1861 order = 1
1862 IF (my_ref) THEN
1863 ! calculate integrals <a|x^n|p>
1864 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true., refpoint=rf, &
1865 particle_set=particle_set, cell=cell)
1866 ELSE
1867 CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.true.)
1868 END IF
1869
1870 ! *** Set up a sorting index
1871 CALL sap_sort(sap_int)
1872
1873 ALLOCATE (basis_set(nkind))
1874 DO ikind = 1, nkind
1875 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1876 IF (ASSOCIATED(orb_basis_set)) THEN
1877 basis_set(ikind)%gto_basis_set => orb_basis_set
1878 ELSE
1879 NULLIFY (basis_set(ikind)%gto_basis_set)
1880 END IF
1881 END DO
1882
1883 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_all)
1884 ! *** All integrals needed have been calculated and stored in sap_int
1885 ! *** We now calculate the commutator matrix elements
1886
1887!$OMP PARALLEL &
1888!$OMP DEFAULT (NONE) &
1889!$OMP SHARED (basis_set, matrix_rv, &
1890!$OMP sap_int, nkind, eps_ppnl, locks, sab_all, &
1891!$OMP particle_set, direction_Or) &
1892!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1893!$OMP iab, irow, icol, blocks_rv, &
1894!$OMP found, iac, ibc, alist_ac, alist_bc, &
1895!$OMP na, np, nb, kkind, kac, kbc, i, lock_num, &
1896!$OMP hash, natom, delta, gamma, achint, bchint, acint, bcint)
1897
1898!$OMP SINGLE
1899!$ ALLOCATE (locks(nlock))
1900!$OMP END SINGLE
1901
1902!$OMP DO
1903!$ DO lock_num = 1, nlock
1904!$ call omp_init_lock(locks(lock_num))
1905!$ END DO
1906!$OMP END DO
1907
1908!$OMP DO SCHEDULE(GUIDED)
1909
1910 DO slot = 1, sab_all(1)%nl_size
1911
1912 ikind = sab_all(1)%nlist_task(slot)%ikind
1913 jkind = sab_all(1)%nlist_task(slot)%jkind
1914 iatom = sab_all(1)%nlist_task(slot)%iatom
1915 jatom = sab_all(1)%nlist_task(slot)%jatom
1916 cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1917 rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1918
1919 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) cycle
1920 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) cycle
1921 iab = ikind + nkind*(jkind - 1)
1922
1923 irow = iatom
1924 icol = jatom
1925
1926 ! allocate blocks
1927 ALLOCATE (blocks_rv(3, 3))
1928
1929 ! get blocks
1930 DO i = 1, 3
1931 DO j = 1, 3
1932 CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, &
1933 blocks_rv(i, j)%block, found)
1934 blocks_rv(i, j)%block(:, :) = 0.0_dp
1935 cpassert(found)
1936 END DO
1937 END DO
1938
1939 ! loop over all kinds for projector atom
1940 ! < iatom | katom > h < katom | jatom >
1941 DO kkind = 1, nkind
1942 iac = ikind + nkind*(kkind - 1)
1943 ibc = jkind + nkind*(kkind - 1)
1944 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
1945 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) cycle
1946 CALL get_alist(sap_int(iac), alist_ac, iatom)
1947 CALL get_alist(sap_int(ibc), alist_bc, jatom)
1948 IF (.NOT. ASSOCIATED(alist_ac)) cycle
1949 IF (.NOT. ASSOCIATED(alist_bc)) cycle
1950 DO kac = 1, alist_ac%nclist
1951 DO kbc = 1, alist_bc%nclist
1952 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
1953
1954 IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1955 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) cycle
1956 acint => alist_ac%clist(kac)%acint
1957 bcint => alist_bc%clist(kbc)%acint
1958 achint => alist_ac%clist(kac)%achint
1959 bchint => alist_bc%clist(kbc)%achint
1960 na = SIZE(acint, 1)
1961 np = SIZE(acint, 2)
1962 nb = SIZE(bcint, 1)
1963!$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1964!$ CALL omp_set_lock(locks(hash))
1965
1966 !! The atom index is alist_ac%clist(kac)%catom
1967 ! The coordinate is particle_set(alist_ac%clist(kac)%catom)%r(:)
1968 IF (direction_or) THEN ! V * r_delta * (R^eta_gamma - R^nu_gamma)
1969 DO delta = 1, 3
1970 DO gamma = 1, 3
1971 blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1972 = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1973 matmul(achint(1:na, 1:np, i_1), transpose(bcint(1:nb, 1:np, delta + 1))) &
1974 *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1975 END DO
1976 END DO
1977 ELSE ! r_delta * V * (R^eta_gamma - R^nu_gamma)
1978 DO delta = 1, 3
1979 DO gamma = 1, 3
1980 blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1981 = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1982 matmul(achint(1:na, 1:np, delta + 1), transpose(bcint(1:nb, 1:np, i_1))) &
1983 *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1984 END DO
1985 END DO
1986 END IF
1987
1988!$ CALL omp_unset_lock(locks(hash))
1989 EXIT ! We have found a match and there can be only one single match
1990 END IF
1991 END DO
1992 END DO
1993 END DO
1994 DO delta = 1, 3
1995 DO gamma = 1, 3
1996 NULLIFY (blocks_rv(gamma, delta)%block)
1997 END DO
1998 END DO
1999 DEALLOCATE (blocks_rv)
2000 END DO
2001
2002!$OMP DO
2003!$ DO lock_num = 1, nlock
2004!$ call omp_destroy_lock(locks(lock_num))
2005!$ END DO
2006!$OMP END DO
2007
2008!$OMP SINGLE
2009!$ DEALLOCATE (locks)
2010!$OMP END SINGLE NOWAIT
2011
2012!$OMP END PARALLEL
2013
2014 CALL release_sap_int(sap_int)
2015
2016 DEALLOCATE (basis_set)
2017
2018 CALL timestop(handle)
2019
2020 END SUBROUTINE build_com_vnl_giao
2021
2022END 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:225
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:74
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:1800
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:1560
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:693
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, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
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, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
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:996
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:550
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:60
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
external_potential_types::gth_potential_p_type
Definition external_potential_types.F:229
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
external_potential_types::sgp_potential_p_type
Definition external_potential_types.F:237
external_potential_types::sgp_potential_type
Definition external_potential_types.F:174
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_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