Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian
9 : !> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
10 : !> \par History
11 : !> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01]
12 : !> - full rewrite [jhu, 2009-01-23]
13 : ! **************************************************************************************************
14 : MODULE commutator_rpnl
15 : USE 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 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief ...
69 : !> \param matrix_rv ...
70 : !> \param qs_kind_set ...
71 : !> \param sab_orb ...
72 : !> \param sap_ppnl ...
73 : !> \param eps_ppnl ...
74 : ! **************************************************************************************************
75 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
91 0 : nsgf_seta
92 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
93 : LOGICAL :: found, gpot, ppnl_present, spot
94 : REAL(KIND=dp) :: dac, ppnl_radius
95 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_ppnl, set_radius_a
99 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, rpgfa, sphi_a, vprj_ppnl, x_block, &
100 0 : y_block, z_block, zeta
101 0 : 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 0 : TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
105 : TYPE(gth_potential_type), POINTER :: gth_potential
106 0 : 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 0 : DIMENSION(:), POINTER :: nl_iterator
110 0 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
111 0 : TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential
112 : TYPE(sgp_potential_type), POINTER :: sgp_potential
113 :
114 0 : CALL timeset(routineN, handle)
115 :
116 0 : ppnl_present = ASSOCIATED(sap_ppnl)
117 :
118 0 : IF (ppnl_present) THEN
119 :
120 0 : 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 0 : maxppnl=maxppnl)
128 :
129 0 : maxl = MAX(maxlgto, maxlppnl)
130 0 : CALL init_orbital_pointers(maxl + 1)
131 :
132 0 : ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl)
133 0 : ldai = ncoset(maxl + 1)
134 :
135 : !sap_int needs to be shared as multiple threads need to access this
136 0 : ALLOCATE (sap_int(nkind*nkind))
137 0 : DO i = 1, nkind*nkind
138 0 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
139 0 : sap_int(i)%nalist = 0
140 : END DO
141 :
142 : !set up direct access to basis and potential
143 0 : ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind))
144 0 : DO ikind = 1, nkind
145 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
146 0 : IF (ASSOCIATED(orb_basis_set)) THEN
147 0 : basis_set(ikind)%gto_basis_set => orb_basis_set
148 : ELSE
149 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
150 : END IF
151 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, &
152 0 : sgp_potential=sgp_potential)
153 0 : IF (ASSOCIATED(gth_potential)) THEN
154 0 : gpotential(ikind)%gth_potential => gth_potential
155 0 : NULLIFY (spotential(ikind)%sgp_potential)
156 0 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
157 0 : spotential(ikind)%sgp_potential => sgp_potential
158 0 : NULLIFY (gpotential(ikind)%gth_potential)
159 : ELSE
160 0 : NULLIFY (gpotential(ikind)%gth_potential)
161 0 : NULLIFY (spotential(ikind)%sgp_potential)
162 : END IF
163 : END DO
164 :
165 0 : maxder = 4
166 : nthread = 1
167 0 : !$ nthread = omp_get_max_threads()
168 :
169 : !calculate the overlap integrals <a|p>
170 0 : 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 0 : !$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 0 : CALL neighbor_list_iterator_release(nl_iterator)
299 :
300 : ! *** Set up a sorting index
301 0 : 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 0 : 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 0 : !$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 0 : CALL neighbor_list_iterator_release(nl_iterator)
402 :
403 0 : CALL release_sap_int(sap_int)
404 :
405 0 : DEALLOCATE (basis_set, gpotential, spotential)
406 :
407 : END IF !ppnl_present
408 :
409 0 : CALL timestop(handle)
410 :
411 0 : 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 196 : SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, &
438 294 : 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 98 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
471 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
472 98 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_rrv, blocks_rrv_vrr, blocks_rv, &
473 98 : blocks_rvr, blocks_rxrv
474 98 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_r_doublecom, blocks_r_rxvr, &
475 98 : blocks_rxvr_r
476 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
477 98 : DIMENSION(:) :: basis_set
478 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
479 98 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
480 :
481 : !$ INTEGER(kind=omp_lock_kind), &
482 98 : !$ ALLOCATABLE, DIMENSION(:) :: locks
483 : !$ INTEGER :: lock_num, hash
484 : !$ INTEGER, PARAMETER :: nlock = 501
485 :
486 98 : ppnl_present = ASSOCIATED(sap_ppnl)
487 98 : IF (.NOT. ppnl_present) RETURN
488 :
489 76 : CALL timeset(routineN, handle)
490 :
491 76 : my_r_doublecom = .FALSE.
492 76 : my_r_rxvr = .FALSE.
493 76 : my_rxvr_r = .FALSE.
494 76 : my_rxrv = .FALSE.
495 76 : my_rrv = .FALSE.
496 76 : my_rv = .FALSE.
497 76 : my_rvr = .FALSE.
498 76 : my_rrv_vrr = .FALSE.
499 76 : IF (PRESENT(matrix_r_doublecom)) my_r_doublecom = .TRUE.
500 76 : IF (PRESENT(matrix_r_rxvr)) my_r_rxvr = .TRUE.
501 76 : IF (PRESENT(matrix_rxvr_r)) my_rxvr_r = .TRUE.
502 76 : IF (PRESENT(matrix_rxrv)) my_rxrv = .TRUE.
503 76 : IF (PRESENT(matrix_rrv)) my_rrv = .TRUE.
504 76 : IF (PRESENT(matrix_rv)) my_rv = .TRUE.
505 76 : IF (PRESENT(matrix_rvr)) my_rvr = .TRUE.
506 76 : IF (PRESENT(matrix_rrv_vrr)) my_rrv_vrr = .TRUE.
507 76 : 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 0 : CPABORT('No dbcsr matrix provided for commutator calculation!')
509 : END IF
510 :
511 76 : natom = SIZE(particle_set)
512 :
513 76 : IF (my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) THEN
514 48 : order = 2
515 48 : CPASSERT(PRESENT(ref_point)) ! need reference point for r x [r,Vnl] and [rr,Vnl]
516 28 : ELSE IF (my_rvr .OR. my_rrv_vrr) THEN
517 6 : order = 2
518 : ELSE
519 22 : order = 1
520 : END IF
521 :
522 : ! When we want the double commutator [[Vnl, r], r], we also want to fix the pseudoatom
523 76 : IF (my_r_doublecom) THEN
524 6 : CPASSERT(PRESENT(pseudoatom))
525 : END IF
526 :
527 184 : periodic = ANY(cell%perd > 0)
528 76 : my_ref = .FALSE.
529 76 : IF (PRESENT(ref_point)) THEN
530 54 : IF (.NOT. periodic) THEN
531 18 : rf = ref_point
532 18 : my_ref = .TRUE.
533 : ELSE ! use my_ref = False in periodic case, corresponds to distributed ref point
534 36 : IF (order > 1) THEN
535 36 : 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 76 : nkind = SIZE(qs_kind_set)
541 :
542 : !sap_int needs to be shared as multiple threads need to access this
543 76 : NULLIFY (sap_int)
544 724 : ALLOCATE (sap_int(nkind*nkind))
545 572 : DO i = 1, nkind*nkind
546 496 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
547 572 : sap_int(i)%nalist = 0
548 : END DO
549 :
550 76 : 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 18 : particle_set=particle_set, cell=cell)
554 : ELSE
555 58 : 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 76 : CALL sap_sort(sap_int)
560 :
561 412 : ALLOCATE (basis_set(nkind))
562 260 : DO ikind = 1, nkind
563 184 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
564 260 : IF (ASSOCIATED(orb_basis_set)) THEN
565 184 : basis_set(ikind)%gto_basis_set => orb_basis_set
566 : ELSE
567 0 : 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 76 : 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 76 : !$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 76 : CALL release_sap_int(sap_int)
1541 :
1542 76 : DEALLOCATE (basis_set)
1543 :
1544 76 : CALL timestop(handle)
1545 :
1546 250 : 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 8 : 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 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1583 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1584 8 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1585 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1586 8 : DIMENSION(:) :: basis_set
1587 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1588 8 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1589 :
1590 : !$ INTEGER(kind=omp_lock_kind), &
1591 8 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1592 : !$ INTEGER :: lock_num, hash
1593 : !$ INTEGER, PARAMETER :: nlock = 501
1594 :
1595 8 : ppnl_present = ASSOCIATED(sap_ppnl)
1596 8 : IF (.NOT. ppnl_present) RETURN
1597 :
1598 8 : CALL timeset(routineN, handle)
1599 :
1600 8 : my_ref = .FALSE.
1601 8 : IF (PRESENT(refpoint)) THEN
1602 8 : my_ref = .TRUE.
1603 8 : CPASSERT(PRESENT(cell))
1604 : END IF
1605 :
1606 8 : natom = SIZE(particle_set)
1607 8 : nkind = SIZE(qs_kind_set)
1608 :
1609 : ! allocate integral storage
1610 8 : NULLIFY (sap_int)
1611 56 : ALLOCATE (sap_int(nkind*nkind))
1612 40 : DO ind = 1, nkind*nkind
1613 32 : NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1614 40 : sap_int(ind)%nalist = 0
1615 : END DO
1616 :
1617 : ! build integrals over GTO + projector functions, refpoint actually
1618 8 : 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 8 : IF (my_ref) THEN
1621 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE., refpoint=refpoint, &
1622 8 : particle_set=particle_set, cell=cell)
1623 : ELSE
1624 0 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE.)
1625 : END IF
1626 :
1627 8 : CALL sap_sort(sap_int)
1628 :
1629 : ! get access to basis sets
1630 40 : ALLOCATE (basis_set(nkind))
1631 24 : DO ikind = 1, nkind
1632 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1633 24 : IF (ASSOCIATED(orb_basis_set)) THEN
1634 16 : basis_set(ikind)%gto_basis_set => orb_basis_set
1635 : ELSE
1636 0 : 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 8 : !$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 8 : DEALLOCATE (basis_set)
1779 8 : CALL release_sap_int(sap_int)
1780 :
1781 8 : CALL timestop(handle)
1782 :
1783 16 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1825 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1826 0 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1827 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1828 0 : DIMENSION(:) :: basis_set
1829 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1830 0 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1831 :
1832 : !$ INTEGER(kind=omp_lock_kind), &
1833 0 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1834 : !$ INTEGER :: lock_num, hash
1835 : !$ INTEGER, PARAMETER :: nlock = 501
1836 :
1837 0 : ppnl_present = ASSOCIATED(sap_ppnl)
1838 0 : IF (.NOT. ppnl_present) RETURN
1839 :
1840 0 : CALL timeset(routineN, handle)
1841 :
1842 0 : natom = SIZE(particle_set)
1843 :
1844 0 : my_ref = .FALSE.
1845 0 : IF (PRESENT(ref_point)) THEN
1846 0 : CPASSERT(PRESENT(cell)) ! need cell as well if refpoint is provided
1847 0 : rf = ref_point
1848 0 : my_ref = .TRUE.
1849 : END IF
1850 :
1851 0 : nkind = SIZE(qs_kind_set)
1852 :
1853 : ! sap_int needs to be shared as multiple threads need to access this
1854 0 : NULLIFY (sap_int)
1855 0 : ALLOCATE (sap_int(nkind*nkind))
1856 0 : DO i = 1, nkind*nkind
1857 0 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1858 0 : sap_int(i)%nalist = 0
1859 : END DO
1860 :
1861 0 : order = 1
1862 0 : 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 0 : particle_set=particle_set, cell=cell)
1866 : ELSE
1867 0 : 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 0 : CALL sap_sort(sap_int)
1872 :
1873 0 : ALLOCATE (basis_set(nkind))
1874 0 : DO ikind = 1, nkind
1875 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1876 0 : IF (ASSOCIATED(orb_basis_set)) THEN
1877 0 : basis_set(ikind)%gto_basis_set => orb_basis_set
1878 : ELSE
1879 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1880 : END IF
1881 : END DO
1882 :
1883 0 : 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 0 : !$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 0 : CALL release_sap_int(sap_int)
2015 :
2016 0 : DEALLOCATE (basis_set)
2017 :
2018 0 : CALL timestop(handle)
2019 :
2020 0 : END SUBROUTINE build_com_vnl_giao
2021 :
2022 : END MODULE commutator_rpnl
|
