Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian
9 : !> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
10 : !> \par History
11 : !> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01]
12 : !> - full rewrite [jhu, 2009-01-23]
13 : ! **************************************************************************************************
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 : ELSEIF (spot) THEN
221 : CPABORT('SGP not implemented')
222 : ELSE
223 : CPABORT('PPNL unknown')
224 : END IF
225 : !$OMP CRITICAL(sap_int_critical)
226 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
227 : sap_int(iac)%a_kind = ikind
228 : sap_int(iac)%p_kind = kkind
229 : sap_int(iac)%nalist = nlist
230 : ALLOCATE (sap_int(iac)%alist(nlist))
231 : DO i = 1, nlist
232 : NULLIFY (sap_int(iac)%alist(i)%clist)
233 : sap_int(iac)%alist(i)%aatom = 0
234 : sap_int(iac)%alist(i)%nclist = 0
235 : END DO
236 : END IF
237 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
238 : sap_int(iac)%alist(ilist)%aatom = iatom
239 : sap_int(iac)%alist(ilist)%nclist = nneighbor
240 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
241 : DO i = 1, nneighbor
242 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
243 : END DO
244 : END IF
245 : !$OMP END CRITICAL(sap_int_critical)
246 : dac = SQRT(SUM(rac*rac))
247 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
248 : clist%catom = katom
249 : clist%cell = cell_c
250 : clist%rac = rac
251 : ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
252 : clist%achint(nsgfa, nppnl, maxder))
253 : clist%acint = 0._dp
254 : clist%achint = 0._dp
255 : clist%nsgf_cnt = 0
256 : NULLIFY (clist%sgf_list)
257 : DO iset = 1, nseta
258 : ncoa = npgfa(iset)*ncoset(la_max(iset))
259 : sgfa = first_sgfa(1, iset)
260 : work = 0._dp
261 : prjc = 1
262 : DO l = 0, lppnl
263 : nprjc = nprj_ppnl(l)*nco(l)
264 : IF (nprjc == 0) CYCLE
265 : rprjc(1) = ppnl_radius
266 : IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE
267 : lc_max = l + 2*(nprj_ppnl(l) - 1)
268 : lc_min = l
269 : zetc(1) = alpha_ppnl(l)
270 : ncoc = ncoset(lc_max)
271 : ! Calculate the primitive overlap and dipole moment integrals
272 : CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
273 : lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab(:, :, 1), 0, .FALSE., ai_work, ldai)
274 : CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
275 : lc_max, 1, zetc, rprjc, 1, rac, (/0._dp, 0._dp, 0._dp/), sab(:, :, 2:4))
276 : ! *** Transformation step projector functions (cartesian->spherical) ***
277 : DO i = 1, maxder
278 : CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, 1, i), ldsab, &
279 : cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, 1, i), ldsab)
280 : END DO
281 : prjc = prjc + nprjc
282 : END DO
283 : DO i = 1, maxder
284 : ! Contraction step (basis functions)
285 : CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
286 : work(1, 1, i), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
287 : ! Multiply with interaction matrix(h)
288 : CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, &
289 : vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa)
290 : END DO
291 : END DO
292 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
293 : clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
294 : END DO
295 :
296 : DEALLOCATE (sab, ai_work, work)
297 : !$OMP END PARALLEL
298 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 116 : SUBROUTINE build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, &
438 174 : 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 58 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
471 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
472 58 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_rrv, blocks_rrv_vrr, blocks_rv, &
473 58 : blocks_rvr, blocks_rxrv
474 58 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_r_doublecom, blocks_r_rxvr, &
475 58 : blocks_rxvr_r
476 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
477 58 : DIMENSION(:) :: basis_set
478 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
479 58 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
480 :
481 : !$ INTEGER(kind=omp_lock_kind), &
482 58 : !$ ALLOCATABLE, DIMENSION(:) :: locks
483 : !$ INTEGER :: lock_num, hash
484 : !$ INTEGER, PARAMETER :: nlock = 501
485 :
486 58 : ppnl_present = ASSOCIATED(sap_ppnl)
487 58 : IF (.NOT. ppnl_present) RETURN
488 :
489 36 : CALL timeset(routineN, handle)
490 :
491 36 : my_r_doublecom = .FALSE.
492 36 : my_r_rxvr = .FALSE.
493 36 : my_rxvr_r = .FALSE.
494 36 : my_rxrv = .FALSE.
495 36 : my_rrv = .FALSE.
496 36 : my_rv = .FALSE.
497 36 : my_rvr = .FALSE.
498 36 : my_rrv_vrr = .FALSE.
499 36 : IF (PRESENT(matrix_r_doublecom)) my_r_doublecom = .TRUE.
500 36 : IF (PRESENT(matrix_r_rxvr)) my_r_rxvr = .TRUE.
501 36 : IF (PRESENT(matrix_rxvr_r)) my_rxvr_r = .TRUE.
502 36 : IF (PRESENT(matrix_rxrv)) my_rxrv = .TRUE.
503 36 : IF (PRESENT(matrix_rrv)) my_rrv = .TRUE.
504 36 : IF (PRESENT(matrix_rv)) my_rv = .TRUE.
505 36 : IF (PRESENT(matrix_rvr)) my_rvr = .TRUE.
506 36 : IF (PRESENT(matrix_rrv_vrr)) my_rrv_vrr = .TRUE.
507 36 : 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 36 : natom = SIZE(particle_set)
512 :
513 36 : IF (my_rxrv .OR. my_rrv .OR. my_r_rxvr .OR. my_rxvr_r .OR. my_r_doublecom) THEN
514 16 : order = 2
515 16 : CPASSERT(PRESENT(ref_point)) ! need reference point for r x [r,Vnl] and [rr,Vnl]
516 20 : ELSE IF (my_rvr .OR. my_rrv_vrr) THEN
517 2 : order = 2
518 : ELSE
519 18 : order = 1
520 : END IF
521 :
522 : ! When we want the double commutator [[Vnl, r], r], we also want to fix the pseudoatom
523 36 : IF (my_r_doublecom) THEN
524 6 : CPASSERT(PRESENT(pseudoatom))
525 : END IF
526 :
527 144 : periodic = ANY(cell%perd > 0)
528 36 : my_ref = .FALSE.
529 36 : IF (PRESENT(ref_point)) THEN
530 18 : 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 0 : IF (order .GT. 1) THEN
535 0 : 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 36 : nkind = SIZE(qs_kind_set)
541 :
542 : !sap_int needs to be shared as multiple threads need to access this
543 36 : NULLIFY (sap_int)
544 252 : ALLOCATE (sap_int(nkind*nkind))
545 180 : DO i = 1, nkind*nkind
546 144 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
547 180 : sap_int(i)%nalist = 0
548 : END DO
549 :
550 36 : 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 18 : 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 36 : CALL sap_sort(sap_int)
560 :
561 180 : ALLOCATE (basis_set(nkind))
562 108 : DO ikind = 1, nkind
563 72 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
564 108 : IF (ASSOCIATED(orb_basis_set)) THEN
565 72 : 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 36 : 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 36 : !$OMP acint, achint, bcint, bchint)
593 :
594 : !$OMP SINGLE
595 : !$ ALLOCATE (locks(nlock))
596 : !$OMP END SINGLE
597 :
598 : !$OMP DO
599 : !$ DO lock_num = 1, nlock
600 : !$ call omp_init_lock(locks(lock_num))
601 : !$ END DO
602 : !$OMP END DO
603 :
604 : !$OMP DO SCHEDULE(GUIDED)
605 :
606 : DO slot = 1, sab_all(1)%nl_size
607 :
608 : ikind = sab_all(1)%nlist_task(slot)%ikind
609 : jkind = sab_all(1)%nlist_task(slot)%jkind
610 : iatom = sab_all(1)%nlist_task(slot)%iatom
611 : jatom = sab_all(1)%nlist_task(slot)%jatom
612 : cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
613 : rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
614 :
615 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
616 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
617 : iab = ikind + nkind*(jkind - 1)
618 :
619 : IF (do_symmetric) THEN
620 : IF (iatom <= jatom) THEN
621 : irow = iatom
622 : icol = jatom
623 : ELSE
624 : irow = jatom
625 : icol = iatom
626 : END IF
627 : ELSE
628 : irow = iatom
629 : icol = jatom
630 : END IF
631 :
632 : ! allocate blocks
633 : IF (my_rv) THEN
634 : ALLOCATE (blocks_rv(3))
635 : END IF
636 : IF (my_rxrv) THEN
637 : ALLOCATE (blocks_rxrv(3))
638 : END IF
639 : IF (my_rrv) THEN
640 : ALLOCATE (blocks_rrv(6))
641 : END IF
642 : IF (my_rvr) THEN
643 : ALLOCATE (blocks_rvr(6))
644 : END IF
645 : IF (my_rrv_vrr) THEN
646 : ALLOCATE (blocks_rrv_vrr(6))
647 : END IF
648 : IF (my_r_rxvr) THEN
649 : ALLOCATE (blocks_r_rxvr(3, 3))
650 : END IF
651 :
652 : IF (my_rxvr_r) THEN
653 : ALLOCATE (blocks_rxvr_r(3, 3))
654 : END IF
655 :
656 : IF (my_r_doublecom) THEN
657 : ALLOCATE (blocks_r_doublecom(3, 3))
658 : END IF
659 :
660 : ! get blocks
661 : IF (my_rv) THEN
662 : DO ind = 1, 3
663 : CALL dbcsr_get_block_p(matrix_rv(ind)%matrix, irow, icol, blocks_rv(ind)%block, found)
664 : END DO
665 : END IF
666 :
667 : IF (my_rxrv) THEN
668 : DO ind = 1, 3
669 : CALL dbcsr_get_block_p(matrix_rxrv(ind)%matrix, irow, icol, blocks_rxrv(ind)%block, found)
670 : blocks_rxrv(ind)%block(:, :) = 0._dp
671 : END DO
672 : END IF
673 :
674 : IF (my_rrv) THEN
675 : DO ind = 1, 6
676 : CALL dbcsr_get_block_p(matrix_rrv(ind)%matrix, irow, icol, blocks_rrv(ind)%block, found)
677 : END DO
678 : END IF
679 :
680 : IF (my_rvr) THEN
681 : DO ind = 1, 6
682 : CALL dbcsr_get_block_p(matrix_rvr(ind)%matrix, irow, icol, blocks_rvr(ind)%block, found)
683 : END DO
684 : END IF
685 :
686 : IF (my_rrv_vrr) THEN
687 : DO ind = 1, 6
688 : CALL dbcsr_get_block_p(matrix_rrv_vrr(ind)%matrix, irow, icol, blocks_rrv_vrr(ind)%block, found)
689 : END DO
690 : END IF
691 :
692 : IF (my_r_rxvr) THEN
693 : DO ind = 1, 3
694 : DO ind2 = 1, 3
695 : CALL dbcsr_get_block_p(matrix_r_rxvr(ind, ind2)%matrix, irow, icol, &
696 : blocks_r_rxvr(ind, ind2)%block, found)
697 : blocks_r_rxvr(ind, ind2)%block(:, :) = 0._dp
698 : END DO
699 : END DO
700 : END IF
701 :
702 : IF (my_rxvr_r) THEN
703 : DO ind = 1, 3
704 : DO ind2 = 1, 3
705 : CALL dbcsr_get_block_p(matrix_rxvr_r(ind, ind2)%matrix, irow, icol, &
706 : blocks_rxvr_r(ind, ind2)%block, found)
707 : blocks_rxvr_r(ind, ind2)%block(:, :) = 0._dp
708 : END DO
709 : END DO
710 : END IF
711 :
712 : IF (my_r_doublecom) THEN
713 : DO ind = 1, 3
714 : DO ind2 = 1, 3
715 : CALL dbcsr_get_block_p(matrix_r_doublecom(ind, ind2)%matrix, irow, icol, &
716 : blocks_r_doublecom(ind, ind2)%block, found)
717 : blocks_r_doublecom(ind, ind2)%block(:, :) = 0._dp
718 : END DO
719 : END DO
720 : END IF
721 :
722 : ! check whether all blocks are associated
723 : go = .TRUE.
724 : IF (my_rv) THEN
725 : asso_rv = (ASSOCIATED(blocks_rv(1)%block) .AND. ASSOCIATED(blocks_rv(2)%block) .AND. &
726 : ASSOCIATED(blocks_rv(3)%block))
727 : go = go .AND. asso_rv
728 : END IF
729 :
730 : IF (my_rxrv) THEN
731 : asso_rxrv = (ASSOCIATED(blocks_rxrv(1)%block) .AND. ASSOCIATED(blocks_rxrv(2)%block) .AND. &
732 : ASSOCIATED(blocks_rxrv(3)%block))
733 : go = go .AND. asso_rxrv
734 : END IF
735 :
736 : IF (my_rrv) THEN
737 : asso_rrv = (ASSOCIATED(blocks_rrv(1)%block) .AND. ASSOCIATED(blocks_rrv(2)%block) .AND. &
738 : ASSOCIATED(blocks_rrv(3)%block) .AND. ASSOCIATED(blocks_rrv(4)%block) .AND. &
739 : ASSOCIATED(blocks_rrv(5)%block) .AND. ASSOCIATED(blocks_rrv(6)%block))
740 : go = go .AND. asso_rrv
741 : END IF
742 :
743 : IF (my_rvr) THEN
744 : asso_rvr = (ASSOCIATED(blocks_rvr(1)%block) .AND. ASSOCIATED(blocks_rvr(2)%block) .AND. &
745 : ASSOCIATED(blocks_rvr(3)%block) .AND. ASSOCIATED(blocks_rvr(4)%block) .AND. &
746 : ASSOCIATED(blocks_rvr(5)%block) .AND. ASSOCIATED(blocks_rvr(6)%block))
747 : go = go .AND. asso_rvr
748 : END IF
749 :
750 : IF (my_rrv_vrr) THEN
751 : asso_rrv_vrr = (ASSOCIATED(blocks_rrv_vrr(1)%block) .AND. ASSOCIATED(blocks_rrv_vrr(2)%block) .AND. &
752 : ASSOCIATED(blocks_rrv_vrr(3)%block) .AND. ASSOCIATED(blocks_rrv_vrr(4)%block) .AND. &
753 : ASSOCIATED(blocks_rrv_vrr(5)%block) .AND. ASSOCIATED(blocks_rrv_vrr(6)%block))
754 : go = go .AND. asso_rrv_vrr
755 : END IF
756 :
757 : IF (my_r_rxvr) THEN
758 : asso_r_rxvr = .TRUE.
759 : DO ind = 1, 3
760 : DO ind2 = 1, 3
761 : asso_r_rxvr = asso_r_rxvr .AND. ASSOCIATED(blocks_r_rxvr(ind, ind2)%block)
762 : END DO
763 : END DO
764 : go = go .AND. asso_r_rxvr
765 : END IF
766 :
767 : IF (my_rxvr_r) THEN
768 : asso_rxvr_r = .TRUE.
769 : DO ind = 1, 3
770 : DO ind2 = 1, 3
771 : asso_rxvr_r = asso_rxvr_r .AND. ASSOCIATED(blocks_rxvr_r(ind, ind2)%block)
772 : END DO
773 : END DO
774 : go = go .AND. asso_rxvr_r
775 : END IF
776 :
777 : IF (my_r_doublecom) THEN
778 : asso_r_doublecom = .TRUE.
779 : DO ind = 1, 3
780 : DO ind2 = 1, 3
781 : asso_r_doublecom = asso_r_doublecom .AND. ASSOCIATED(blocks_r_doublecom(ind, ind2)%block)
782 : END DO
783 : END DO
784 : go = go .AND. asso_r_doublecom
785 : END IF
786 :
787 : ! loop over all kinds for projector atom
788 : ! < iatom | katom > h < katom | jatom >
789 : IF (go) THEN
790 : DO kkind = 1, nkind
791 : iac = ikind + nkind*(kkind - 1)
792 : ibc = jkind + nkind*(kkind - 1)
793 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
794 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
795 : CALL get_alist(sap_int(iac), alist_ac, iatom)
796 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
797 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
798 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
799 : DO kac = 1, alist_ac%nclist
800 : DO kbc = 1, alist_bc%nclist
801 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
802 : IF (PRESENT(pseudoatom)) THEN
803 : IF (alist_ac%clist(kac)%catom /= pseudoatom) CYCLE
804 : END IF
805 :
806 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
807 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
808 : acint => alist_ac%clist(kac)%acint
809 : bcint => alist_bc%clist(kbc)%acint
810 : achint => alist_ac%clist(kac)%achint
811 : bchint => alist_bc%clist(kbc)%achint
812 : na = SIZE(acint, 1)
813 : np = SIZE(acint, 2)
814 : nb = SIZE(bcint, 1)
815 : !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
816 : !$ CALL omp_set_lock(locks(hash))
817 : IF (my_rv) THEN
818 : ! r*Vnl
819 : ! with LAPACK
820 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 2), na, &
821 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(1)%block, SIZE(blocks_rv(1)%block, 1)) ! xV
822 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 3), na, &
823 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(2)%block, SIZE(blocks_rv(2)%block, 1)) ! yV
824 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 4), na, &
825 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rv(3)%block, SIZE(blocks_rv(3)%block, 1)) ! zV
826 : IF (iatom <= jatom) THEN
827 : ! with MATMUL
828 : blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) + &
829 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xV
830 : blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) + &
831 : MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yV
832 : blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) + &
833 : MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zV
834 : ELSE
835 : blocks_rv(1)%block(1:nb, 1:na) = blocks_rv(1)%block(1:nb, 1:na) + &
836 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))
837 : blocks_rv(2)%block(1:nb, 1:na) = blocks_rv(2)%block(1:nb, 1:na) + &
838 : MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))
839 : blocks_rv(3)%block(1:nb, 1:na) = blocks_rv(3)%block(1:nb, 1:na) + &
840 : MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 1)))
841 : END IF
842 : ! -Vnl r
843 : ! with LAPACK
844 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
845 : ! bcint(1, 1, 2), nb, 1.0_dp, blocks_rv(1)%block, SIZE(blocks_rv(1)%block, 1)) ! -Vx
846 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
847 : ! bcint(1, 1, 3), nb, 1.0_dp, blocks_rv(2)%block, SIZE(blocks_rv(2)%block, 1)) ! -Vy
848 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
849 : ! bcint(1, 1, 4), nb, 1.0_dp, blocks_rv(3)%block, SIZE(blocks_rv(3)%block, 1)) ! -Vz
850 : ! with MATMUL
851 : IF (iatom <= jatom) THEN
852 : blocks_rv(1)%block(1:na, 1:nb) = blocks_rv(1)%block(1:na, 1:nb) - &
853 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! -Vx
854 : blocks_rv(2)%block(1:na, 1:nb) = blocks_rv(2)%block(1:na, 1:nb) - &
855 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! -Vy
856 : blocks_rv(3)%block(1:na, 1:nb) = blocks_rv(3)%block(1:na, 1:nb) - &
857 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! -Vz
858 : ELSE
859 : blocks_rv(1)%block(1:nb, 1:na) = blocks_rv(1)%block(1:nb, 1:na) - &
860 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 2)))
861 : blocks_rv(2)%block(1:nb, 1:na) = blocks_rv(2)%block(1:nb, 1:na) - &
862 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 3)))
863 : blocks_rv(3)%block(1:nb, 1:na) = blocks_rv(3)%block(1:nb, 1:na) - &
864 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 4)))
865 : END IF
866 :
867 : END IF
868 :
869 : IF (my_rxrv) THEN
870 : ! x-component (y [z,Vnl] - z [y, Vnl])
871 : ! with LAPACK
872 : ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 9), na, &
873 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! yzV
874 : ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 3), na, &
875 : ! bcint(1, 1, 4), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! -yVz
876 : ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 9), na, &
877 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! -zyV
878 : ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 4), na, &
879 : ! bcint(1, 1, 3), nb, 1.0_dp, blocks_rxrv(1)%block, SIZE(blocks_rxrv(1)%block, 1)) ! zVy
880 : ! with MATMUL
881 : IF (iatom <= jatom) THEN
882 : blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
883 : MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yzV
884 : blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
885 : MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! -yVz
886 : blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) - &
887 : MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! -zyV
888 : blocks_rxrv(1)%block(1:na, 1:nb) = blocks_rxrv(1)%block(1:na, 1:nb) + &
889 : MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! zVy
890 : ELSE
891 : blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
892 : MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! yzV
893 : blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
894 : MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 4))) ! -yVz
895 : blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) - &
896 : MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! -zyV
897 : blocks_rxrv(1)%block(1:nb, 1:na) = blocks_rxrv(1)%block(1:nb, 1:na) + &
898 : MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 3))) ! zVy
899 : END IF
900 :
901 : ! y-component (z [x,Vnl] - x [z, Vnl])
902 : ! with LAPACK
903 : ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 7), na, &
904 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! zxV
905 : ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 4), na, &
906 : ! bcint(1, 1, 2), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! -zVx
907 : ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 7), na, &
908 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! -xzV
909 : ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 2), na, &
910 : ! bcint(1, 1, 4), nb, 1.0_dp, blocks_rxrv(2)%block, SIZE(blocks_rxrv(2)%block, 1)) ! xVz
911 : ! with MATMUL
912 : IF (iatom <= jatom) THEN
913 : blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
914 : MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zxV
915 : blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
916 : MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! -zVx
917 : blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) - &
918 : MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! -xzV
919 : blocks_rxrv(2)%block(1:na, 1:nb) = blocks_rxrv(2)%block(1:na, 1:nb) + &
920 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! xVz
921 : ELSE
922 : blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
923 : MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! zxV
924 : blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
925 : MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 2))) ! -zVx
926 : blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) - &
927 : MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! -xzV
928 : blocks_rxrv(2)%block(1:nb, 1:na) = blocks_rxrv(2)%block(1:nb, 1:na) + &
929 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 4))) ! xVz
930 : END IF
931 :
932 : ! z-component (x [y,Vnl] - y [x, Vnl])
933 : ! with LAPACK
934 : ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 6), na, &
935 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! xyV
936 : ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 2), na, &
937 : ! bcint(1, 1, 3), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! -xVy
938 : ! CALL dgemm("N", "T", na, nb, np, -1.0_dp, achint(1, 1, 6), na, &
939 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! -yxV
940 : ! CALL dgemm("N", "T", na, nb, np, 1.0_dp, achint(1, 1, 3), na, &
941 : ! bcint(1, 1, 2), nb, 1.0_dp, blocks_rxrv(3)%block, SIZE(blocks_rxrv(3)%block, 1)) ! yVx
942 : ! with MATMUL
943 : IF (iatom <= jatom) THEN
944 : blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
945 : MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xyV
946 : blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
947 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! -xVy
948 : blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) - &
949 : MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! -yxV
950 : blocks_rxrv(3)%block(1:na, 1:nb) = blocks_rxrv(3)%block(1:na, 1:nb) + &
951 : MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! zVx
952 : ELSE
953 : blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
954 : MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! xyV
955 : blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
956 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 3))) ! -xVy
957 : blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) - &
958 : MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! -yxV
959 : blocks_rxrv(3)%block(1:nb, 1:na) = blocks_rxrv(3)%block(1:nb, 1:na) + &
960 : MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 2))) ! zVx
961 : END IF
962 : END IF
963 :
964 : IF (my_rrv) THEN
965 : ! r_alpha * r_beta * Vnl
966 : ! with LAPACK
967 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 5), na, &
968 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(1)%block, SIZE(blocks_rrv(1)%block, 1)) ! xxV
969 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 6), na, &
970 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(2)%block, SIZE(blocks_rrv(2)%block, 1)) ! xyV
971 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 7), na, &
972 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(3)%block, SIZE(blocks_rrv(3)%block, 1)) ! xzV
973 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 8), na, &
974 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(4)%block, SIZE(blocks_rrv(4)%block, 1)) ! yyV
975 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 9), na, &
976 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(5)%block, SIZE(blocks_rrv(5)%block, 1)) ! yzV
977 : ! CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 10), na, &
978 : ! bcint(1, 1, 1), nb, 1.0_dp, blocks_rrv(6)%block, SIZE(blocks_rrv(6)%block, 1)) ! zzV
979 : ! with MATMUL
980 : IF (iatom <= jatom) THEN
981 : blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) + &
982 : MATMUL(achint(1:na, 1:np, 5), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xxV
983 : blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) + &
984 : MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xyV
985 : blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) + &
986 : MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xzV
987 : blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) + &
988 : MATMUL(achint(1:na, 1:np, 8), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yyV
989 : blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) + &
990 : MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yzV
991 : blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) + &
992 : MATMUL(achint(1:na, 1:np, 10), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zzV
993 : ELSE
994 : blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) + &
995 : MATMUL(bchint(1:nb, 1:np, 5), TRANSPOSE(acint(1:na, 1:np, 1))) ! xxV
996 : blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) + &
997 : MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! xyV
998 : blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) + &
999 : MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! xzV
1000 : blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) + &
1001 : MATMUL(bchint(1:nb, 1:np, 8), TRANSPOSE(acint(1:na, 1:np, 1))) ! yyV
1002 : blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) + &
1003 : MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! yzV
1004 : blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) + &
1005 : MATMUL(bchint(1:nb, 1:np, 10), TRANSPOSE(acint(1:na, 1:np, 1))) ! zzV
1006 : END IF
1007 :
1008 : ! - Vnl * r_alpha * r_beta
1009 : ! with LAPACK
1010 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1011 : ! bcint(1, 1, 5), nb, 1.0_dp, blocks_rrv(1)%block, SIZE(blocks_rrv(1)%block, 1)) ! Vxx
1012 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1013 : ! bcint(1, 1, 6), nb, 1.0_dp, blocks_rrv(2)%block, SIZE(blocks_rrv(2)%block, 1)) ! Vxy
1014 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1015 : ! bcint(1, 1, 7), nb, 1.0_dp, blocks_rrv(3)%block, SIZE(blocks_rrv(3)%block, 1)) ! Vxz
1016 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1017 : ! bcint(1, 1, 8), nb, 1.0_dp, blocks_rrv(4)%block, SIZE(blocks_rrv(4)%block, 1)) ! Vyy
1018 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1019 : ! bcint(1, 1, 9), nb, 1.0_dp, blocks_rrv(5)%block, SIZE(blocks_rrv(5)%block, 1)) ! Vyz
1020 : ! CALL dgemm("N", "T", na, nb, np, -1._dp, achint(1, 1, 1), na, &
1021 : ! bcint(1, 1, 10), nb, 1.0_dp, blocks_rrv(6)%block, SIZE(blocks_rrv(6)%block, 1)) ! Vzz
1022 : ! with MATMUL
1023 : IF (iatom <= jatom) THEN
1024 : blocks_rrv(1)%block(1:na, 1:nb) = blocks_rrv(1)%block(1:na, 1:nb) - &
1025 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 5))) ! -Vxx
1026 : blocks_rrv(2)%block(1:na, 1:nb) = blocks_rrv(2)%block(1:na, 1:nb) - &
1027 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 6))) ! -Vxy
1028 : blocks_rrv(3)%block(1:na, 1:nb) = blocks_rrv(3)%block(1:na, 1:nb) - &
1029 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 7))) ! -Vxz
1030 : blocks_rrv(4)%block(1:na, 1:nb) = blocks_rrv(4)%block(1:na, 1:nb) - &
1031 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 8))) ! -Vyy
1032 : blocks_rrv(5)%block(1:na, 1:nb) = blocks_rrv(5)%block(1:na, 1:nb) - &
1033 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 9))) ! -Vyz
1034 : blocks_rrv(6)%block(1:na, 1:nb) = blocks_rrv(6)%block(1:na, 1:nb) - &
1035 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 10))) ! -Vzz
1036 : ELSE
1037 : blocks_rrv(1)%block(1:nb, 1:na) = blocks_rrv(1)%block(1:nb, 1:na) - &
1038 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 5))) ! -Vxx
1039 : blocks_rrv(2)%block(1:nb, 1:na) = blocks_rrv(2)%block(1:nb, 1:na) - &
1040 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 6))) ! -Vxy
1041 : blocks_rrv(3)%block(1:nb, 1:na) = blocks_rrv(3)%block(1:nb, 1:na) - &
1042 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 7))) ! -Vxz
1043 : blocks_rrv(4)%block(1:nb, 1:na) = blocks_rrv(4)%block(1:nb, 1:na) - &
1044 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 8))) ! -Vyy
1045 : blocks_rrv(5)%block(1:nb, 1:na) = blocks_rrv(5)%block(1:nb, 1:na) - &
1046 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 9))) ! -Vyz
1047 : blocks_rrv(6)%block(1:nb, 1:na) = blocks_rrv(6)%block(1:nb, 1:na) - &
1048 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 10))) ! -Vzz
1049 : END IF
1050 : END IF
1051 :
1052 : IF (my_rvr) THEN
1053 : ! r_alpha * Vnl * r_beta
1054 : IF (iatom <= jatom) THEN
1055 : blocks_rvr(1)%block(1:na, 1:nb) = blocks_rvr(1)%block(1:na, 1:nb) + &
1056 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 2))) ! xVx
1057 : blocks_rvr(2)%block(1:na, 1:nb) = blocks_rvr(2)%block(1:na, 1:nb) + &
1058 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! xVy
1059 : blocks_rvr(3)%block(1:na, 1:nb) = blocks_rvr(3)%block(1:na, 1:nb) + &
1060 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! xVz
1061 : blocks_rvr(4)%block(1:na, 1:nb) = blocks_rvr(4)%block(1:na, 1:nb) + &
1062 : MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 3))) ! yVy
1063 : blocks_rvr(5)%block(1:na, 1:nb) = blocks_rvr(5)%block(1:na, 1:nb) + &
1064 : MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! yVz
1065 : blocks_rvr(6)%block(1:na, 1:nb) = blocks_rvr(6)%block(1:na, 1:nb) + &
1066 : MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 4))) ! zVz
1067 : ELSE
1068 : blocks_rvr(1)%block(1:nb, 1:na) = blocks_rvr(1)%block(1:nb, 1:na) + &
1069 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 2))) ! xVx
1070 : blocks_rvr(2)%block(1:nb, 1:na) = blocks_rvr(2)%block(1:nb, 1:na) + &
1071 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 3))) ! xVy
1072 : blocks_rvr(3)%block(1:nb, 1:na) = blocks_rvr(3)%block(1:nb, 1:na) + &
1073 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 4))) ! xVz
1074 : blocks_rvr(4)%block(1:nb, 1:na) = blocks_rvr(4)%block(1:nb, 1:na) + &
1075 : MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 3))) ! yVy
1076 : blocks_rvr(5)%block(1:nb, 1:na) = blocks_rvr(5)%block(1:nb, 1:na) + &
1077 : MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 4))) ! yVz
1078 : blocks_rvr(6)%block(1:nb, 1:na) = blocks_rvr(6)%block(1:nb, 1:na) + &
1079 : MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 4))) ! zVz
1080 : END IF
1081 : END IF
1082 :
1083 : IF (my_rrv_vrr) THEN
1084 : ! r_alpha * r_beta * Vnl
1085 : IF (iatom <= jatom) THEN
1086 : blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
1087 : MATMUL(achint(1:na, 1:np, 5), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xxV
1088 : blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
1089 : MATMUL(achint(1:na, 1:np, 6), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xyV
1090 : blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
1091 : MATMUL(achint(1:na, 1:np, 7), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! xzV
1092 : blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
1093 : MATMUL(achint(1:na, 1:np, 8), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yyV
1094 : blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
1095 : MATMUL(achint(1:na, 1:np, 9), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! yzV
1096 : blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
1097 : MATMUL(achint(1:na, 1:np, 10), TRANSPOSE(bcint(1:nb, 1:np, 1))) ! zzV
1098 : ELSE
1099 : blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
1100 : MATMUL(bchint(1:nb, 1:np, 5), TRANSPOSE(acint(1:na, 1:np, 1))) ! xxV
1101 : blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
1102 : MATMUL(bchint(1:nb, 1:np, 6), TRANSPOSE(acint(1:na, 1:np, 1))) ! xyV
1103 : blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
1104 : MATMUL(bchint(1:nb, 1:np, 7), TRANSPOSE(acint(1:na, 1:np, 1))) ! xzV
1105 : blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
1106 : MATMUL(bchint(1:nb, 1:np, 8), TRANSPOSE(acint(1:na, 1:np, 1))) ! yyV
1107 : blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
1108 : MATMUL(bchint(1:nb, 1:np, 9), TRANSPOSE(acint(1:na, 1:np, 1))) ! yzV
1109 : blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
1110 : MATMUL(bchint(1:nb, 1:np, 10), TRANSPOSE(acint(1:na, 1:np, 1))) ! zzV
1111 : END IF
1112 : ! + Vnl * r_alpha * r_beta
1113 : IF (iatom <= jatom) THEN
1114 : blocks_rrv_vrr(1)%block(1:na, 1:nb) = blocks_rrv_vrr(1)%block(1:na, 1:nb) + &
1115 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 5))) ! +Vxx
1116 : blocks_rrv_vrr(2)%block(1:na, 1:nb) = blocks_rrv_vrr(2)%block(1:na, 1:nb) + &
1117 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 6))) ! +Vxy
1118 : blocks_rrv_vrr(3)%block(1:na, 1:nb) = blocks_rrv_vrr(3)%block(1:na, 1:nb) + &
1119 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 7))) ! +Vxz
1120 : blocks_rrv_vrr(4)%block(1:na, 1:nb) = blocks_rrv_vrr(4)%block(1:na, 1:nb) + &
1121 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 8))) ! +Vyy
1122 : blocks_rrv_vrr(5)%block(1:na, 1:nb) = blocks_rrv_vrr(5)%block(1:na, 1:nb) + &
1123 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 9))) ! +Vyz
1124 : blocks_rrv_vrr(6)%block(1:na, 1:nb) = blocks_rrv_vrr(6)%block(1:na, 1:nb) + &
1125 : MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 10))) ! +Vzz
1126 : ELSE
1127 : blocks_rrv_vrr(1)%block(1:nb, 1:na) = blocks_rrv_vrr(1)%block(1:nb, 1:na) + &
1128 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 5))) ! +Vxx
1129 : blocks_rrv_vrr(2)%block(1:nb, 1:na) = blocks_rrv_vrr(2)%block(1:nb, 1:na) + &
1130 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 6))) ! +Vxy
1131 : blocks_rrv_vrr(3)%block(1:nb, 1:na) = blocks_rrv_vrr(3)%block(1:nb, 1:na) + &
1132 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 7))) ! +Vxz
1133 : blocks_rrv_vrr(4)%block(1:nb, 1:na) = blocks_rrv_vrr(4)%block(1:nb, 1:na) + &
1134 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 8))) ! +Vyy
1135 : blocks_rrv_vrr(5)%block(1:nb, 1:na) = blocks_rrv_vrr(5)%block(1:nb, 1:na) + &
1136 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 9))) ! +Vyz
1137 : blocks_rrv_vrr(6)%block(1:nb, 1:na) = blocks_rrv_vrr(6)%block(1:nb, 1:na) + &
1138 : MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 10))) ! +Vzz
1139 : END IF
1140 : END IF
1141 :
1142 : ! The indices are stored in i_1, i_x, ..., i_zzz
1143 : ! matrix_r_rxvr(alpha, beta)
1144 : ! = sum_(gamma delta) epsilon_(alpha gamma delta)
1145 : ! (r_beta * r_gamma * V_nl * r_delta - r_beta * r_gamma * r_delta * V_nl)
1146 : ! = sum_(gamma delta) epsilon_(alpha gamma delta) r_beta * r_gamma * V_nl * r_delta
1147 :
1148 : ! TODO: is this set to zero before?
1149 : IF (my_r_rxvr) THEN
1150 : ! beta = 1
1151 : ! matrix_r_rxvr(x, x) = x * y * V_nl * z - x * z * V_nl * y
1152 : blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
1153 : blocks_r_rxvr(1, 1)%block(1:na, 1:nb) + &
1154 : MATMUL(achint(1:na, 1:np, i_xy), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1155 : blocks_r_rxvr(1, 1)%block(1:na, 1:nb) = &
1156 : blocks_r_rxvr(1, 1)%block(1:na, 1:nb) - &
1157 : MATMUL(achint(1:na, 1:np, i_xz), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1158 :
1159 : ! matrix_r_rxvr(y, x) = x * z * V_nl * x - x * x * V_nl * z
1160 : blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
1161 : blocks_r_rxvr(2, 1)%block(1:na, 1:nb) + &
1162 : MATMUL(achint(1:na, 1:np, i_xz), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1163 : blocks_r_rxvr(2, 1)%block(1:na, 1:nb) = &
1164 : blocks_r_rxvr(2, 1)%block(1:na, 1:nb) - &
1165 : MATMUL(achint(1:na, 1:np, i_xx), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1166 :
1167 : ! matrix_r_rxvr(z, x) = x * x * V_nl * y - x * y * V_nl * x
1168 : blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
1169 : blocks_r_rxvr(3, 1)%block(1:na, 1:nb) + &
1170 : MATMUL(achint(1:na, 1:np, i_xx), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1171 : blocks_r_rxvr(3, 1)%block(1:na, 1:nb) = &
1172 : blocks_r_rxvr(3, 1)%block(1:na, 1:nb) - &
1173 : MATMUL(achint(1:na, 1:np, i_xy), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1174 :
1175 : ! beta = 2
1176 : ! matrix_r_rxvr(x, y) = y * y * V_nl * z - y * z * V_nl * y
1177 : blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
1178 : blocks_r_rxvr(1, 2)%block(1:na, 1:nb) + &
1179 : MATMUL(achint(1:na, 1:np, i_yy), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1180 : blocks_r_rxvr(1, 2)%block(1:na, 1:nb) = &
1181 : blocks_r_rxvr(1, 2)%block(1:na, 1:nb) - &
1182 : MATMUL(achint(1:na, 1:np, i_yz), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1183 :
1184 : ! matrix_r_rxvr(y, y) = y * z * V_nl * x - y * x * V_nl * z
1185 : blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
1186 : blocks_r_rxvr(2, 2)%block(1:na, 1:nb) + &
1187 : MATMUL(achint(1:na, 1:np, i_yz), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1188 : blocks_r_rxvr(2, 2)%block(1:na, 1:nb) = &
1189 : blocks_r_rxvr(2, 2)%block(1:na, 1:nb) - &
1190 : MATMUL(achint(1:na, 1:np, i_yx), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1191 :
1192 : ! matrix_r_rxvr(z, y) = y * x * V_nl * y - y * y * V_nl * x
1193 : blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
1194 : blocks_r_rxvr(3, 2)%block(1:na, 1:nb) + &
1195 : MATMUL(achint(1:na, 1:np, i_yx), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1196 : blocks_r_rxvr(3, 2)%block(1:na, 1:nb) = &
1197 : blocks_r_rxvr(3, 2)%block(1:na, 1:nb) - &
1198 : MATMUL(achint(1:na, 1:np, i_yy), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1199 :
1200 : ! beta = 3
1201 : ! matrix_r_rxvr(x, z) = z * y * V_nl * z - z * z * V_nl * y
1202 : blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
1203 : blocks_r_rxvr(1, 3)%block(1:na, 1:nb) + &
1204 : MATMUL(achint(1:na, 1:np, i_zy), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1205 : blocks_r_rxvr(1, 3)%block(1:na, 1:nb) = &
1206 : blocks_r_rxvr(1, 3)%block(1:na, 1:nb) - &
1207 : MATMUL(achint(1:na, 1:np, i_zz), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1208 :
1209 : ! matrix_r_rxvr(y, z) = z * z * V_nl * x - z * x * V_nl * z
1210 : blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
1211 : blocks_r_rxvr(2, 3)%block(1:na, 1:nb) + &
1212 : MATMUL(achint(1:na, 1:np, i_zz), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1213 : blocks_r_rxvr(2, 3)%block(1:na, 1:nb) = &
1214 : blocks_r_rxvr(2, 3)%block(1:na, 1:nb) - &
1215 : MATMUL(achint(1:na, 1:np, i_zx), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1216 :
1217 : ! matrix_r_rxvr(z, z) = z * x * V_nl * y - z * y * V_nl * x
1218 : blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
1219 : blocks_r_rxvr(3, 3)%block(1:na, 1:nb) + &
1220 : MATMUL(achint(1:na, 1:np, i_zx), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1221 : blocks_r_rxvr(3, 3)%block(1:na, 1:nb) = &
1222 : blocks_r_rxvr(3, 3)%block(1:na, 1:nb) - &
1223 : MATMUL(achint(1:na, 1:np, i_zy), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1224 :
1225 : END IF ! my_r_rxvr
1226 :
1227 : ! The indices are stored in i_1, i_x, ..., i_zzz
1228 : ! This will put into blocks_rxvr_r
1229 : ! matrix_rxvr_r(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
1230 : ! r_gamma * V_nl * r_delta * r_beta
1231 : IF (my_rxvr_r) THEN
1232 : ! beta = 1
1233 : ! matrix_rxvr_r(x, x) = yV zx - zV yx
1234 : blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
1235 : blocks_rxvr_r(1, 1)%block(1:na, 1:nb) + &
1236 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_zx)))
1237 : blocks_rxvr_r(1, 1)%block(1:na, 1:nb) = &
1238 : blocks_rxvr_r(1, 1)%block(1:na, 1:nb) - &
1239 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_yx)))
1240 :
1241 : ! matrix_rxvr_r(y, x) = zV xx - xV zx
1242 : blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
1243 : blocks_rxvr_r(2, 1)%block(1:na, 1:nb) + &
1244 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_xx)))
1245 : blocks_rxvr_r(2, 1)%block(1:na, 1:nb) = &
1246 : blocks_rxvr_r(2, 1)%block(1:na, 1:nb) - &
1247 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_zx)))
1248 :
1249 : ! matrix_rxvr_r(z, x) = xV yx - yV xx
1250 : blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
1251 : blocks_rxvr_r(3, 1)%block(1:na, 1:nb) + &
1252 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_yx)))
1253 : blocks_rxvr_r(3, 1)%block(1:na, 1:nb) = &
1254 : blocks_rxvr_r(3, 1)%block(1:na, 1:nb) - &
1255 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_xx)))
1256 :
1257 : ! beta = 2
1258 : ! matrix_rxvr_r(x, y) = yV zy - zV yy
1259 : blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
1260 : blocks_rxvr_r(1, 2)%block(1:na, 1:nb) + &
1261 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_zy)))
1262 : blocks_rxvr_r(1, 2)%block(1:na, 1:nb) = &
1263 : blocks_rxvr_r(1, 2)%block(1:na, 1:nb) - &
1264 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_yy)))
1265 :
1266 : ! matrix_rxvr_r(y, y) = zV xy - xV zy
1267 : blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
1268 : blocks_rxvr_r(2, 2)%block(1:na, 1:nb) + &
1269 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_xy)))
1270 : blocks_rxvr_r(2, 2)%block(1:na, 1:nb) = &
1271 : blocks_rxvr_r(2, 2)%block(1:na, 1:nb) - &
1272 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_zy)))
1273 :
1274 : ! matrix_rxvr_r(z, y) = xV yy - yV xy
1275 : blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
1276 : blocks_rxvr_r(3, 2)%block(1:na, 1:nb) + &
1277 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_yy)))
1278 : blocks_rxvr_r(3, 2)%block(1:na, 1:nb) = &
1279 : blocks_rxvr_r(3, 2)%block(1:na, 1:nb) - &
1280 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_xy)))
1281 :
1282 : ! beta = 3
1283 : ! matrix_rxvr_r(x, z) = yV zz - zV yz
1284 : blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
1285 : blocks_rxvr_r(1, 3)%block(1:na, 1:nb) + &
1286 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_zz)))
1287 : blocks_rxvr_r(1, 3)%block(1:na, 1:nb) = &
1288 : blocks_rxvr_r(1, 3)%block(1:na, 1:nb) - &
1289 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_yz)))
1290 :
1291 : ! matrix_rxvr_r(y, z) = zV xz - xV zz
1292 : blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
1293 : blocks_rxvr_r(2, 3)%block(1:na, 1:nb) + &
1294 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_xz)))
1295 : blocks_rxvr_r(2, 3)%block(1:na, 1:nb) = &
1296 : blocks_rxvr_r(2, 3)%block(1:na, 1:nb) - &
1297 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_zz)))
1298 :
1299 : ! matrix_rxvr_r(z, z) = xV yz - yV xz
1300 : blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
1301 : blocks_rxvr_r(3, 3)%block(1:na, 1:nb) + &
1302 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_yz)))
1303 : blocks_rxvr_r(3, 3)%block(1:na, 1:nb) = &
1304 : blocks_rxvr_r(3, 3)%block(1:na, 1:nb) - &
1305 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_xz)))
1306 :
1307 : END IF ! my_rxvr_r
1308 :
1309 : ! matrix_r_doublecom(alpha, beta) = sum_(gamma delta) epsilon_(alpha gamma delta)
1310 : ! gamma V^pseudoatom beta delta - gamma beta V^pseudoatom delta
1311 :
1312 : IF (my_r_doublecom) THEN
1313 : ! beta = 1
1314 : ! matrix_r_doublecom(x, x) = yV xz - zV xy - yxV z + zxV y
1315 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1316 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
1317 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_xz)))
1318 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1319 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
1320 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_xy)))
1321 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1322 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) - &
1323 : MATMUL(achint(1:na, 1:np, i_yx), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1324 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) = &
1325 : blocks_r_doublecom(1, 1)%block(1:na, 1:nb) + &
1326 : MATMUL(achint(1:na, 1:np, i_zx), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1327 :
1328 : ! matrix_r_doublecom(y, x) = zV xx - xV xz - zxV x + xxV z
1329 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1330 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
1331 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_xx)))
1332 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1333 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
1334 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_xz)))
1335 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1336 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) - &
1337 : MATMUL(achint(1:na, 1:np, i_zx), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1338 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) = &
1339 : blocks_r_doublecom(2, 1)%block(1:na, 1:nb) + &
1340 : MATMUL(achint(1:na, 1:np, i_xx), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1341 :
1342 : ! matrix_r_doublecom(z, x) = xV xy - yV xx - xxV y + yxV x
1343 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1344 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1345 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_xy)))
1346 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1347 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1348 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_xx)))
1349 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1350 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) - &
1351 : MATMUL(achint(1:na, 1:np, i_xx), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1352 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) = &
1353 : blocks_r_doublecom(3, 1)%block(1:na, 1:nb) + &
1354 : MATMUL(achint(1:na, 1:np, i_yx), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1355 :
1356 : ! beta = 2
1357 : ! matrix_r_doublecom(x, y) = yV yz - zV yy - yyV z + zyV y
1358 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1359 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1360 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_yz)))
1361 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1362 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1363 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_yy)))
1364 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1365 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) - &
1366 : MATMUL(achint(1:na, 1:np, i_yy), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1367 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) = &
1368 : blocks_r_doublecom(1, 2)%block(1:na, 1:nb) + &
1369 : MATMUL(achint(1:na, 1:np, i_zy), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1370 :
1371 : ! matrix_r_doublecom(y, y) = zV yx - xV yz - zyV x + xyV z
1372 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1373 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1374 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_yx)))
1375 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1376 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1377 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_yz)))
1378 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1379 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) - &
1380 : MATMUL(achint(1:na, 1:np, i_zy), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1381 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) = &
1382 : blocks_r_doublecom(2, 2)%block(1:na, 1:nb) + &
1383 : MATMUL(achint(1:na, 1:np, i_xy), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1384 :
1385 : ! matrix_r_doublecom(z, y) = xV yy - yV yx - xyV y + yyV x
1386 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1387 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1388 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_yy)))
1389 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1390 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1391 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_yx)))
1392 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1393 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) - &
1394 : MATMUL(achint(1:na, 1:np, i_xy), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1395 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) = &
1396 : blocks_r_doublecom(3, 2)%block(1:na, 1:nb) + &
1397 : MATMUL(achint(1:na, 1:np, i_yy), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1398 :
1399 : ! beta = 3
1400 : ! matrix_r_doublecom(x, z) = yV zz - zV zy - yzV z + zzV y
1401 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1402 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1403 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_zz)))
1404 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1405 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1406 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_zy)))
1407 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1408 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) - &
1409 : MATMUL(achint(1:na, 1:np, i_yz), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1410 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) = &
1411 : blocks_r_doublecom(1, 3)%block(1:na, 1:nb) + &
1412 : MATMUL(achint(1:na, 1:np, i_zz), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1413 :
1414 : ! matrix_r_doublecom(y, z) = zV zx - xV zz - zzV x + xzV z
1415 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1416 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1417 : MATMUL(achint(1:na, 1:np, i_z), TRANSPOSE(bcint(1:nb, 1:np, i_zx)))
1418 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1419 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1420 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_zz)))
1421 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1422 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) - &
1423 : MATMUL(achint(1:na, 1:np, i_zz), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1424 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) = &
1425 : blocks_r_doublecom(2, 3)%block(1:na, 1:nb) + &
1426 : MATMUL(achint(1:na, 1:np, i_xz), TRANSPOSE(bcint(1:nb, 1:np, i_z)))
1427 :
1428 : ! matrix_r_doublecom(z, z) = xV zy - yV zx - xzV y + yzV x
1429 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1430 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1431 : MATMUL(achint(1:na, 1:np, i_x), TRANSPOSE(bcint(1:nb, 1:np, i_zy)))
1432 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1433 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1434 : MATMUL(achint(1:na, 1:np, i_y), TRANSPOSE(bcint(1:nb, 1:np, i_zx)))
1435 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1436 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) - &
1437 : MATMUL(achint(1:na, 1:np, i_xz), TRANSPOSE(bcint(1:nb, 1:np, i_y)))
1438 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) = &
1439 : blocks_r_doublecom(3, 3)%block(1:na, 1:nb) + &
1440 : MATMUL(achint(1:na, 1:np, i_yz), TRANSPOSE(bcint(1:nb, 1:np, i_x)))
1441 :
1442 : END IF ! my_r_doublecom
1443 : !$ CALL omp_unset_lock(locks(hash))
1444 : EXIT ! We have found a match and there can be only one single match
1445 : END IF
1446 : END DO
1447 : END DO
1448 : END DO
1449 : END IF
1450 : IF (my_rv) THEN
1451 : DO ind = 1, 3
1452 : NULLIFY (blocks_rv(ind)%block)
1453 : END DO
1454 : DEALLOCATE (blocks_rv)
1455 : END IF
1456 : IF (my_rxrv) THEN
1457 : DO ind = 1, 3
1458 : NULLIFY (blocks_rxrv(ind)%block)
1459 : END DO
1460 : DEALLOCATE (blocks_rxrv)
1461 : END IF
1462 : IF (my_rrv) THEN
1463 : DO ind = 1, 6
1464 : NULLIFY (blocks_rrv(ind)%block)
1465 : END DO
1466 : DEALLOCATE (blocks_rrv)
1467 : END IF
1468 : IF (my_rvr) THEN
1469 : DO ind = 1, 6
1470 : NULLIFY (blocks_rvr(ind)%block)
1471 : END DO
1472 : DEALLOCATE (blocks_rvr)
1473 : END IF
1474 : IF (my_rrv_vrr) THEN
1475 : DO ind = 1, 6
1476 : NULLIFY (blocks_rrv_vrr(ind)%block)
1477 : END DO
1478 : DEALLOCATE (blocks_rrv_vrr)
1479 : END IF
1480 : IF (my_r_rxvr) THEN
1481 : DO ind = 1, 3
1482 : DO ind2 = 1, 3
1483 : NULLIFY (blocks_r_rxvr(ind, ind2)%block)
1484 : END DO
1485 : END DO
1486 : DEALLOCATE (blocks_r_rxvr)
1487 : END IF
1488 : IF (my_rxvr_r) THEN
1489 : DO ind = 1, 3
1490 : DO ind2 = 1, 3
1491 : NULLIFY (blocks_rxvr_r(ind, ind2)%block)
1492 : END DO
1493 : END DO
1494 : DEALLOCATE (blocks_rxvr_r)
1495 : END IF
1496 : IF (my_r_doublecom) THEN
1497 : DO ind = 1, 3
1498 : DO ind2 = 1, 3
1499 : NULLIFY (blocks_r_doublecom(ind, ind2)%block)
1500 : END DO
1501 : END DO
1502 : DEALLOCATE (blocks_r_doublecom)
1503 : END IF
1504 : END DO
1505 :
1506 : !$OMP DO
1507 : !$ DO lock_num = 1, nlock
1508 : !$ call omp_destroy_lock(locks(lock_num))
1509 : !$ END DO
1510 : !$OMP END DO
1511 :
1512 : !$OMP SINGLE
1513 : !$ DEALLOCATE (locks)
1514 : !$OMP END SINGLE NOWAIT
1515 :
1516 : !$OMP END PARALLEL
1517 :
1518 36 : CALL release_sap_int(sap_int)
1519 :
1520 36 : DEALLOCATE (basis_set)
1521 :
1522 36 : CALL timestop(handle)
1523 :
1524 130 : END SUBROUTINE build_com_mom_nl
1525 :
1526 : ! **************************************************************************************************
1527 : !> \brief calculate \sum_R_ps (R_ps - R_nu) x [V_nl, r] summing over all pseudized atoms R
1528 : !> \param qs_kind_set ...
1529 : !> \param sab_all ...
1530 : !> \param sap_ppnl ...
1531 : !> \param eps_ppnl ...
1532 : !> \param particle_set ...
1533 : !> \param matrix_mag_nl ...
1534 : !> \param refpoint ...
1535 : !> \param cell ...
1536 : ! **************************************************************************************************
1537 8 : SUBROUTINE build_com_nl_mag(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_mag_nl, refpoint, cell)
1538 :
1539 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1540 : POINTER :: qs_kind_set
1541 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1542 : INTENT(IN), POINTER :: sab_all, sap_ppnl
1543 : REAL(KIND=dp), INTENT(IN) :: eps_ppnl
1544 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1545 : POINTER :: particle_set
1546 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1547 : POINTER :: matrix_mag_nl
1548 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: refpoint
1549 : TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1550 :
1551 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_com_nl_mag'
1552 :
1553 : INTEGER :: handle, iab, iac, iatom, ibc, icol, &
1554 : ikind, ind, irow, jatom, jkind, kac, &
1555 : kbc, kkind, na, natom, nb, nkind, np, &
1556 : order, slot
1557 : INTEGER, DIMENSION(3) :: cell_b
1558 : LOGICAL :: found, go, my_ref, ppnl_present
1559 : REAL(KIND=dp), DIMENSION(3) :: r_b, r_ps, rab
1560 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1561 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1562 8 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: blocks_mag
1563 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1564 8 : DIMENSION(:) :: basis_set
1565 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1566 8 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1567 :
1568 : !$ INTEGER(kind=omp_lock_kind), &
1569 8 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1570 : !$ INTEGER :: lock_num, hash
1571 : !$ INTEGER, PARAMETER :: nlock = 501
1572 :
1573 8 : ppnl_present = ASSOCIATED(sap_ppnl)
1574 8 : IF (.NOT. ppnl_present) RETURN
1575 :
1576 8 : CALL timeset(routineN, handle)
1577 :
1578 8 : my_ref = .FALSE.
1579 8 : IF (PRESENT(refpoint)) THEN
1580 8 : my_ref = .TRUE.
1581 8 : CPASSERT(PRESENT(cell))
1582 : END IF
1583 :
1584 8 : natom = SIZE(particle_set)
1585 8 : nkind = SIZE(qs_kind_set)
1586 :
1587 : ! allocate integral storage
1588 8 : NULLIFY (sap_int)
1589 56 : ALLOCATE (sap_int(nkind*nkind))
1590 40 : DO ind = 1, nkind*nkind
1591 32 : NULLIFY (sap_int(ind)%alist, sap_int(ind)%asort, sap_int(ind)%aindex)
1592 40 : sap_int(ind)%nalist = 0
1593 : END DO
1594 :
1595 : ! build integrals over GTO + projector functions, refpoint actually
1596 8 : order = 1 ! only need first moments (x, y, z)
1597 : ! refpoint actually does not matter in this case, i. e. (order = 1 .and. commutator)
1598 8 : IF (my_ref) THEN
1599 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE., refpoint=refpoint, &
1600 8 : particle_set=particle_set, cell=cell)
1601 : ELSE
1602 0 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE.)
1603 : END IF
1604 :
1605 8 : CALL sap_sort(sap_int)
1606 :
1607 : ! get access to basis sets
1608 40 : ALLOCATE (basis_set(nkind))
1609 24 : DO ikind = 1, nkind
1610 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1611 24 : IF (ASSOCIATED(orb_basis_set)) THEN
1612 16 : basis_set(ikind)%gto_basis_set => orb_basis_set
1613 : ELSE
1614 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1615 : END IF
1616 : END DO
1617 :
1618 : !$OMP PARALLEL &
1619 : !$OMP DEFAULT (NONE) &
1620 : !$OMP SHARED (basis_set, matrix_mag_nl, sap_int, natom, nkind, eps_ppnl, locks, sab_all, &
1621 : !$OMP particle_set, my_ref, refpoint) &
1622 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, lock_num, &
1623 : !$OMP iab, irow, icol, blocks_mag, r_ps, r_b, go, hash, &
1624 : !$OMP found, iac, ibc, alist_ac, alist_bc, acint, bcint, &
1625 8 : !$OMP achint, bchint, na, np, nb, kkind, kac, kbc)
1626 :
1627 : !$OMP SINGLE
1628 : !$ ALLOCATE (locks(nlock))
1629 : !$OMP END SINGLE
1630 :
1631 : !$OMP DO
1632 : !$ DO lock_num = 1, nlock
1633 : !$ call omp_init_lock(locks(lock_num))
1634 : !$ END DO
1635 : !$OMP END DO
1636 :
1637 : !$OMP DO SCHEDULE(GUIDED)
1638 : DO slot = 1, sab_all(1)%nl_size
1639 : ! get indices
1640 : ikind = sab_all(1)%nlist_task(slot)%ikind
1641 : jkind = sab_all(1)%nlist_task(slot)%jkind
1642 : iatom = sab_all(1)%nlist_task(slot)%iatom
1643 : jatom = sab_all(1)%nlist_task(slot)%jatom
1644 : cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1645 : rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1646 :
1647 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
1648 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
1649 : iab = ikind + nkind*(jkind - 1)
1650 :
1651 : IF (iatom <= jatom) THEN
1652 : irow = iatom
1653 : icol = jatom
1654 : ELSE
1655 : irow = jatom
1656 : icol = iatom
1657 : END IF
1658 :
1659 : ! get blocks
1660 : ALLOCATE (blocks_mag(3))
1661 : DO ind = 1, 3
1662 : CALL dbcsr_get_block_p(matrix_mag_nl(ind)%matrix, irow, icol, blocks_mag(ind)%block, found)
1663 : END DO
1664 :
1665 : go = (ASSOCIATED(blocks_mag(1)%block) .AND. ASSOCIATED(blocks_mag(2)%block) .AND. ASSOCIATED(blocks_mag(3)%block))
1666 :
1667 : IF (go) THEN
1668 : DO kkind = 1, nkind
1669 : iac = ikind + nkind*(kkind - 1)
1670 : ibc = jkind + nkind*(kkind - 1)
1671 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
1672 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
1673 : CALL get_alist(sap_int(iac), alist_ac, iatom)
1674 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
1675 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
1676 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
1677 : DO kac = 1, alist_ac%nclist
1678 : DO kbc = 1, alist_bc%nclist
1679 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
1680 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1681 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
1682 :
1683 : acint => alist_ac%clist(kac)%acint
1684 : bcint => alist_bc%clist(kbc)%acint
1685 : achint => alist_ac%clist(kac)%achint
1686 : bchint => alist_bc%clist(kbc)%achint
1687 : na = SIZE(acint, 1)
1688 : np = SIZE(acint, 2)
1689 : nb = SIZE(bcint, 1)
1690 : ! Position of the pseudized atom
1691 : r_ps = particle_set(alist_ac%clist(kac)%catom)%r
1692 : r_b = refpoint
1693 :
1694 : !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1695 : !$ CALL omp_set_lock(locks(hash))
1696 : ! assemble integrals
1697 : IF (iatom <= jatom) THEN
1698 : blocks_mag(1)%block(1:na, 1:nb) = blocks_mag(1)%block(1:na, 1:nb) + &
1699 : (r_ps(2) - r_b(2))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 4))) - &
1700 : MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 1)))) & ! R_y [V_nl, z]
1701 : - (r_ps(3) - r_b(3))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 3))) - &
1702 : MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))) ! - R_z [V_nl, y]
1703 : blocks_mag(2)%block(1:na, 1:nb) = blocks_mag(2)%block(1:na, 1:nb) + &
1704 : (r_ps(3) - r_b(3))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 2))) - &
1705 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))) & ! R_z [V_nl, x]
1706 : - (r_ps(1) - r_b(1))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 4))) - &
1707 : MATMUL(achint(1:na, 1:np, 4), TRANSPOSE(bcint(1:nb, 1:np, 1)))) ! - R_x [V_nl, z]
1708 : blocks_mag(3)%block(1:na, 1:nb) = blocks_mag(3)%block(1:na, 1:nb) + &
1709 : (r_ps(1) - r_b(1))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 3))) - &
1710 : MATMUL(achint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))) & ! R_x [V_nl, y]
1711 : - (r_ps(2) - r_b(2))*(MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 2))) - &
1712 : MATMUL(achint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))) ! - R_y [V_nl, x]
1713 : ELSE
1714 : blocks_mag(1)%block(1:nb, 1:na) = blocks_mag(1)%block(1:nb, 1:na) + &
1715 : (r_ps(2) - r_b(2))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 4))) - &
1716 : MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 1)))) & ! R_y [V_nl, z]
1717 : - (r_ps(3) - r_b(3))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 3))) - &
1718 : MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))) ! - R_z [V_nl, y]
1719 : blocks_mag(2)%block(1:nb, 1:na) = blocks_mag(2)%block(1:nb, 1:na) + &
1720 : (r_ps(3) - r_b(3))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 2))) - &
1721 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))) & ! R_z [V_nl, x]
1722 : - (r_ps(1) - r_b(1))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 4))) - &
1723 : MATMUL(bchint(1:nb, 1:np, 4), TRANSPOSE(acint(1:na, 1:np, 1)))) ! - R_x [V_nl, z]
1724 : blocks_mag(3)%block(1:nb, 1:na) = blocks_mag(3)%block(1:nb, 1:na) + &
1725 : (r_ps(1) - r_b(1))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 3))) - &
1726 : MATMUL(bchint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))) & ! R_x [V_nl, y]
1727 : - (r_ps(2) - r_b(2))*(MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 2))) - &
1728 : MATMUL(bchint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))) ! - R_y [V_nl, x]
1729 : END IF
1730 : !$ CALL omp_unset_lock(locks(hash))
1731 : EXIT ! We have found a match and there can be only one single match
1732 : END IF
1733 : END DO
1734 : END DO
1735 : END DO
1736 : END IF
1737 :
1738 : DO ind = 1, 3
1739 : NULLIFY (blocks_mag(ind)%block)
1740 : END DO
1741 : DEALLOCATE (blocks_mag)
1742 : END DO
1743 :
1744 : !$OMP DO
1745 : !$ DO lock_num = 1, nlock
1746 : !$ call omp_destroy_lock(locks(lock_num))
1747 : !$ END DO
1748 : !$OMP END DO
1749 :
1750 : !$OMP SINGLE
1751 : !$ DEALLOCATE (locks)
1752 : !$OMP END SINGLE NOWAIT
1753 :
1754 : !$OMP END PARALLEL
1755 :
1756 8 : DEALLOCATE (basis_set)
1757 8 : CALL release_sap_int(sap_int)
1758 :
1759 8 : CALL timestop(handle)
1760 :
1761 16 : END SUBROUTINE build_com_nl_mag
1762 :
1763 : ! **************************************************************************************************
1764 : !> \brief Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs
1765 : !> \param qs_kind_set ...
1766 : !> \param sab_all ...
1767 : !> \param sap_ppnl ...
1768 : !> \param eps_ppnl ...
1769 : !> \param particle_set ...
1770 : !> \param matrix_rv ...
1771 : !> \param ref_point ...
1772 : !> \param cell ...
1773 : !> \param direction_Or If set to true: calculate Vnl * r_delta
1774 : !> Otherwise calculate r_delta * Vnl
1775 : ! **************************************************************************************************
1776 0 : SUBROUTINE build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, &
1777 : matrix_rv, ref_point, cell, direction_Or)
1778 :
1779 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1780 : POINTER :: qs_kind_set
1781 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1782 : INTENT(IN), POINTER :: sab_all, sap_ppnl
1783 : REAL(KIND=dp), INTENT(IN) :: eps_ppnl
1784 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1785 : POINTER :: particle_set
1786 : TYPE(dbcsr_p_type), DIMENSION(:, :), &
1787 : INTENT(INOUT), OPTIONAL, POINTER :: matrix_rv
1788 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: ref_point
1789 : TYPE(cell_type), INTENT(IN), OPTIONAL, POINTER :: cell
1790 : LOGICAL :: direction_Or
1791 :
1792 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_com_vnl_giao'
1793 : INTEGER, PARAMETER :: i_1 = 1
1794 :
1795 : INTEGER :: delta, gamma, handle, i, iab, iac, &
1796 : iatom, ibc, icol, ikind, irow, j, &
1797 : jatom, jkind, kac, kbc, kkind, na, &
1798 : natom, nb, nkind, np, order, slot
1799 : INTEGER, DIMENSION(3) :: cell_b
1800 : LOGICAL :: found, my_ref, ppnl_present
1801 : REAL(KIND=dp), DIMENSION(3) :: rab, rf
1802 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint
1803 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
1804 0 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: blocks_rv
1805 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1806 0 : DIMENSION(:) :: basis_set
1807 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1808 0 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
1809 :
1810 : !$ INTEGER(kind=omp_lock_kind), &
1811 0 : !$ ALLOCATABLE, DIMENSION(:) :: locks
1812 : !$ INTEGER :: lock_num, hash
1813 : !$ INTEGER, PARAMETER :: nlock = 501
1814 :
1815 0 : ppnl_present = ASSOCIATED(sap_ppnl)
1816 0 : IF (.NOT. ppnl_present) RETURN
1817 :
1818 0 : CALL timeset(routineN, handle)
1819 :
1820 0 : natom = SIZE(particle_set)
1821 :
1822 0 : my_ref = .FALSE.
1823 0 : IF (PRESENT(ref_point)) THEN
1824 0 : CPASSERT(PRESENT(cell)) ! need cell as well if refpoint is provided
1825 0 : rf = ref_point
1826 0 : my_ref = .TRUE.
1827 : END IF
1828 :
1829 0 : nkind = SIZE(qs_kind_set)
1830 :
1831 : ! sap_int needs to be shared as multiple threads need to access this
1832 0 : NULLIFY (sap_int)
1833 0 : ALLOCATE (sap_int(nkind*nkind))
1834 0 : DO i = 1, nkind*nkind
1835 0 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
1836 0 : sap_int(i)%nalist = 0
1837 : END DO
1838 :
1839 0 : order = 1
1840 0 : IF (my_ref) THEN
1841 : ! calculate integrals <a|x^n|p>
1842 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE., refpoint=rf, &
1843 0 : particle_set=particle_set, cell=cell)
1844 : ELSE
1845 0 : CALL build_sap_ints(sap_int, sap_ppnl, qs_kind_set, order, moment_mode=.TRUE.)
1846 : END IF
1847 :
1848 : ! *** Set up a sorting index
1849 0 : CALL sap_sort(sap_int)
1850 :
1851 0 : ALLOCATE (basis_set(nkind))
1852 0 : DO ikind = 1, nkind
1853 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1854 0 : IF (ASSOCIATED(orb_basis_set)) THEN
1855 0 : basis_set(ikind)%gto_basis_set => orb_basis_set
1856 : ELSE
1857 0 : NULLIFY (basis_set(ikind)%gto_basis_set)
1858 : END IF
1859 : END DO
1860 :
1861 0 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_all)
1862 : ! *** All integrals needed have been calculated and stored in sap_int
1863 : ! *** We now calculate the commutator matrix elements
1864 :
1865 : !$OMP PARALLEL &
1866 : !$OMP DEFAULT (NONE) &
1867 : !$OMP SHARED (basis_set, matrix_rv, &
1868 : !$OMP sap_int, nkind, eps_ppnl, locks, sab_all, &
1869 : !$OMP particle_set, direction_Or) &
1870 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
1871 : !$OMP iab, irow, icol, blocks_rv, &
1872 : !$OMP found, iac, ibc, alist_ac, alist_bc, &
1873 : !$OMP na, np, nb, kkind, kac, kbc, i, lock_num, &
1874 0 : !$OMP hash, natom, delta, gamma, achint, bchint, acint, bcint)
1875 :
1876 : !$OMP SINGLE
1877 : !$ ALLOCATE (locks(nlock))
1878 : !$OMP END SINGLE
1879 :
1880 : !$OMP DO
1881 : !$ DO lock_num = 1, nlock
1882 : !$ call omp_init_lock(locks(lock_num))
1883 : !$ END DO
1884 : !$OMP END DO
1885 :
1886 : !$OMP DO SCHEDULE(GUIDED)
1887 :
1888 : DO slot = 1, sab_all(1)%nl_size
1889 :
1890 : ikind = sab_all(1)%nlist_task(slot)%ikind
1891 : jkind = sab_all(1)%nlist_task(slot)%jkind
1892 : iatom = sab_all(1)%nlist_task(slot)%iatom
1893 : jatom = sab_all(1)%nlist_task(slot)%jatom
1894 : cell_b(:) = sab_all(1)%nlist_task(slot)%cell(:)
1895 : rab(1:3) = sab_all(1)%nlist_task(slot)%r(1:3)
1896 :
1897 : IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE
1898 : IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE
1899 : iab = ikind + nkind*(jkind - 1)
1900 :
1901 : irow = iatom
1902 : icol = jatom
1903 :
1904 : ! allocate blocks
1905 : ALLOCATE (blocks_rv(3, 3))
1906 :
1907 : ! get blocks
1908 : DO i = 1, 3
1909 : DO j = 1, 3
1910 : CALL dbcsr_get_block_p(matrix_rv(i, j)%matrix, irow, icol, &
1911 : blocks_rv(i, j)%block, found)
1912 : blocks_rv(i, j)%block(:, :) = 0.0_dp
1913 : CPASSERT(found)
1914 : END DO
1915 : END DO
1916 :
1917 : ! loop over all kinds for projector atom
1918 : ! < iatom | katom > h < katom | jatom >
1919 : DO kkind = 1, nkind
1920 : iac = ikind + nkind*(kkind - 1)
1921 : ibc = jkind + nkind*(kkind - 1)
1922 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
1923 : IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE
1924 : CALL get_alist(sap_int(iac), alist_ac, iatom)
1925 : CALL get_alist(sap_int(ibc), alist_bc, jatom)
1926 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
1927 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
1928 : DO kac = 1, alist_ac%nclist
1929 : DO kbc = 1, alist_bc%nclist
1930 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
1931 :
1932 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
1933 : IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE
1934 : acint => alist_ac%clist(kac)%acint
1935 : bcint => alist_bc%clist(kbc)%acint
1936 : achint => alist_ac%clist(kac)%achint
1937 : bchint => alist_bc%clist(kbc)%achint
1938 : na = SIZE(acint, 1)
1939 : np = SIZE(acint, 2)
1940 : nb = SIZE(bcint, 1)
1941 : !$ hash = MOD((iatom - 1)*natom + jatom, nlock) + 1
1942 : !$ CALL omp_set_lock(locks(hash))
1943 :
1944 : !! The atom index is alist_ac%clist(kac)%catom
1945 : ! The coordinate is particle_set(alist_ac%clist(kac)%catom)%r(:)
1946 : IF (direction_Or) THEN ! V * r_delta * (R^eta_gamma - R^nu_gamma)
1947 : DO delta = 1, 3
1948 : DO gamma = 1, 3
1949 : blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1950 : = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1951 : MATMUL(achint(1:na, 1:np, i_1), TRANSPOSE(bcint(1:nb, 1:np, delta + 1))) &
1952 : *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1953 : END DO
1954 : END DO
1955 : ELSE ! r_delta * V * (R^eta_gamma - R^nu_gamma)
1956 : DO delta = 1, 3
1957 : DO gamma = 1, 3
1958 : blocks_rv(gamma, delta)%block(1:na, 1:nb) &
1959 : = blocks_rv(gamma, delta)%block(1:na, 1:nb) + &
1960 : MATMUL(achint(1:na, 1:np, delta + 1), TRANSPOSE(bcint(1:nb, 1:np, i_1))) &
1961 : *(particle_set(alist_ac%clist(kac)%catom)%r(gamma) - particle_set(jatom)%r(gamma))
1962 : END DO
1963 : END DO
1964 : END IF
1965 :
1966 : !$ CALL omp_unset_lock(locks(hash))
1967 : EXIT ! We have found a match and there can be only one single match
1968 : END IF
1969 : END DO
1970 : END DO
1971 : END DO
1972 : DO delta = 1, 3
1973 : DO gamma = 1, 3
1974 : NULLIFY (blocks_rv(gamma, delta)%block)
1975 : END DO
1976 : END DO
1977 : DEALLOCATE (blocks_rv)
1978 : END DO
1979 :
1980 : !$OMP DO
1981 : !$ DO lock_num = 1, nlock
1982 : !$ call omp_destroy_lock(locks(lock_num))
1983 : !$ END DO
1984 : !$OMP END DO
1985 :
1986 : !$OMP SINGLE
1987 : !$ DEALLOCATE (locks)
1988 : !$OMP END SINGLE NOWAIT
1989 :
1990 : !$OMP END PARALLEL
1991 :
1992 0 : CALL release_sap_int(sap_int)
1993 :
1994 0 : DEALLOCATE (basis_set)
1995 :
1996 0 : CALL timestop(handle)
1997 :
1998 0 : END SUBROUTINE build_com_vnl_giao
1999 :
2000 : END MODULE commutator_rpnl
|
