Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief Methods used with 3-center overlap type integrals containers
9 : !> \par History
10 : !> - none
11 : !> - 11.2018 fixed OMP race condition in contract3_o3c routine (A.Bussy)
12 : ! **************************************************************************************************
13 : MODULE qs_o3c_methods
14 : USE ai_contraction_sphi, ONLY: abc_contract
15 : USE ai_overlap3, ONLY: overlap3
16 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
17 : gto_basis_set_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
19 : dbcsr_p_type,&
20 : dbcsr_type
21 : USE kinds, ONLY: dp
22 : USE orbital_pointers, ONLY: ncoset
23 : USE qs_o3c_types, ONLY: &
24 : get_o3c_container, get_o3c_iterator_info, get_o3c_vec, o3c_container_type, o3c_iterate, &
25 : o3c_iterator_create, o3c_iterator_release, o3c_iterator_type, o3c_vec_type, &
26 : set_o3c_container
27 :
28 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_o3c_methods'
36 :
37 : PUBLIC :: calculate_o3c_integrals, contract12_o3c, contract3_o3c
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief ...
43 : !> \param o3c ...
44 : !> \param calculate_forces ...
45 : !> \param matrix_p ...
46 : ! **************************************************************************************************
47 0 : SUBROUTINE calculate_o3c_integrals(o3c, calculate_forces, matrix_p)
48 : TYPE(o3c_container_type), POINTER :: o3c
49 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
50 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
51 : POINTER :: matrix_p
52 :
53 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_o3c_integrals'
54 :
55 : INTEGER :: egfa, egfb, egfc, handle, i, iatom, icol, ikind, irow, iset, ispin, j, jatom, &
56 : jkind, jset, katom, kkind, kset, mepos, ncoa, ncob, ncoc, ni, nj, nk, nseta, nsetb, &
57 : nsetc, nspin, nthread, sgfa, sgfb, sgfc
58 0 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
59 0 : lc_min, npgfa, npgfb, npgfc, nsgfa, &
60 0 : nsgfb, nsgfc
61 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfc
62 : LOGICAL :: do_force, found, trans
63 : REAL(KIND=dp) :: dij, dik, djk, fpre
64 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pmat
65 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sabc, sabc_contr
66 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: iabdc, iadbc, idabc, sabdc, sdabc
67 : REAL(KIND=dp), DIMENSION(3) :: rij, rik, rjk
68 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_c
69 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fi, fj, fk, pblock, rpgfa, rpgfb, rpgfc, &
70 0 : sphi_a, sphi_b, sphi_c, tvec, zeta, &
71 0 : zetb, zetc
72 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: iabc
73 0 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b, &
74 0 : basis_set_list_c
75 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, basis_set_c
76 : TYPE(o3c_iterator_type) :: o3c_iterator
77 :
78 0 : CALL timeset(routineN, handle)
79 :
80 0 : do_force = .FALSE.
81 0 : IF (PRESENT(calculate_forces)) do_force = calculate_forces
82 0 : CALL get_o3c_container(o3c, nspin=nspin)
83 :
84 : ! basis sets
85 : CALL get_o3c_container(o3c, basis_set_list_a=basis_set_list_a, &
86 0 : basis_set_list_b=basis_set_list_b, basis_set_list_c=basis_set_list_c)
87 :
88 : nthread = 1
89 0 : !$ nthread = omp_get_max_threads()
90 0 : CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
91 :
92 : !$OMP PARALLEL DEFAULT(NONE) &
93 : !$OMP SHARED (nthread,o3c_iterator,ncoset,nspin,basis_set_list_a,basis_set_list_b,&
94 : !$OMP basis_set_list_c,do_force,matrix_p)&
95 : !$OMP PRIVATE (mepos,ikind,jkind,kkind,basis_set_a,basis_set_b,basis_set_c,rij,rik,rjk,&
96 : !$OMP first_sgfa,la_max,la_min,npgfa,nseta,nsgfa,rpgfa,set_radius_a,sphi_a,zeta,&
97 : !$OMP first_sgfb,lb_max,lb_min,npgfb,nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,&
98 : !$OMP first_sgfc,lc_max,lc_min,npgfc,nsetc,nsgfc,rpgfc,set_radius_c,sphi_c,zetc,&
99 : !$OMP iset,jset,kset,dij,dik,djk,ni,nj,nk,iabc,idabc,iadbc,iabdc,tvec,fi,fj,fk,ncoa,&
100 : !$OMP ncob,ncoc,sabc,sabc_contr,sdabc,sabdc,sgfa,sgfb,sgfc,egfa,egfb,egfc,i,j,&
101 0 : !$OMP pblock,pmat,ispin,iatom,jatom,katom,irow,icol,found,trans,fpre)
102 :
103 : mepos = 0
104 : !$ mepos = omp_get_thread_num()
105 :
106 : DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
107 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, &
108 : ikind=ikind, jkind=jkind, kkind=kkind, rij=rij, rik=rik, &
109 : integral=iabc, tvec=tvec, force_i=fi, force_j=fj, force_k=fk)
110 : CPASSERT(.NOT. ASSOCIATED(iabc))
111 : CPASSERT(.NOT. ASSOCIATED(tvec))
112 : CPASSERT(.NOT. ASSOCIATED(fi))
113 : CPASSERT(.NOT. ASSOCIATED(fj))
114 : CPASSERT(.NOT. ASSOCIATED(fk))
115 : ! basis
116 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
117 : basis_set_b => basis_set_list_b(jkind)%gto_basis_set
118 : basis_set_c => basis_set_list_c(kkind)%gto_basis_set
119 : ! center A
120 : first_sgfa => basis_set_a%first_sgf
121 : la_max => basis_set_a%lmax
122 : la_min => basis_set_a%lmin
123 : npgfa => basis_set_a%npgf
124 : nseta = basis_set_a%nset
125 : nsgfa => basis_set_a%nsgf_set
126 : rpgfa => basis_set_a%pgf_radius
127 : set_radius_a => basis_set_a%set_radius
128 : sphi_a => basis_set_a%sphi
129 : zeta => basis_set_a%zet
130 : ! center B
131 : first_sgfb => basis_set_b%first_sgf
132 : lb_max => basis_set_b%lmax
133 : lb_min => basis_set_b%lmin
134 : npgfb => basis_set_b%npgf
135 : nsetb = basis_set_b%nset
136 : nsgfb => basis_set_b%nsgf_set
137 : rpgfb => basis_set_b%pgf_radius
138 : set_radius_b => basis_set_b%set_radius
139 : sphi_b => basis_set_b%sphi
140 : zetb => basis_set_b%zet
141 : ! center C (RI)
142 : first_sgfc => basis_set_c%first_sgf
143 : lc_max => basis_set_c%lmax
144 : lc_min => basis_set_c%lmin
145 : npgfc => basis_set_c%npgf
146 : nsetc = basis_set_c%nset
147 : nsgfc => basis_set_c%nsgf_set
148 : rpgfc => basis_set_c%pgf_radius
149 : set_radius_c => basis_set_c%set_radius
150 : sphi_c => basis_set_c%sphi
151 : zetc => basis_set_c%zet
152 :
153 : ni = SUM(nsgfa)
154 : nj = SUM(nsgfb)
155 : nk = SUM(nsgfc)
156 :
157 : ALLOCATE (iabc(ni, nj, nk))
158 : iabc(:, :, :) = 0.0_dp
159 : IF (do_force) THEN
160 : ALLOCATE (fi(nk, 3), fj(nk, 3), fk(nk, 3))
161 : fi(:, :) = 0.0_dp
162 : fj(:, :) = 0.0_dp
163 : fk(:, :) = 0.0_dp
164 : ALLOCATE (idabc(ni, nj, nk, 3))
165 : idabc(:, :, :, :) = 0.0_dp
166 : ALLOCATE (iadbc(ni, nj, nk, 3))
167 : iadbc(:, :, :, :) = 0.0_dp
168 : ALLOCATE (iabdc(ni, nj, nk, 3))
169 : iabdc(:, :, :, :) = 0.0_dp
170 : ELSE
171 : NULLIFY (fi, fj, fk)
172 : END IF
173 : ALLOCATE (tvec(nk, nspin))
174 : tvec(:, :) = 0.0_dp
175 :
176 : rjk(1:3) = rik(1:3) - rij(1:3)
177 : dij = NORM2(rij)
178 : dik = NORM2(rik)
179 : djk = NORM2(rjk)
180 :
181 : DO iset = 1, nseta
182 : DO jset = 1, nsetb
183 : IF (set_radius_a(iset) + set_radius_b(jset) < dij) CYCLE
184 : DO kset = 1, nsetc
185 : IF (set_radius_a(iset) + set_radius_c(kset) < dik) CYCLE
186 : IF (set_radius_b(jset) + set_radius_c(kset) < djk) CYCLE
187 :
188 : ncoa = npgfa(iset)*ncoset(la_max(iset))
189 : ncob = npgfb(jset)*ncoset(lb_max(jset))
190 : ncoc = npgfc(kset)*ncoset(lc_max(kset))
191 :
192 : sgfa = first_sgfa(1, iset)
193 : sgfb = first_sgfb(1, jset)
194 : sgfc = first_sgfc(1, kset)
195 :
196 : egfa = sgfa + nsgfa(iset) - 1
197 : egfb = sgfb + nsgfb(jset) - 1
198 : egfc = sgfc + nsgfc(kset) - 1
199 :
200 : IF (ncoa*ncob*ncoc > 0) THEN
201 : ALLOCATE (sabc(ncoa, ncob, ncoc))
202 : sabc(:, :, :) = 0.0_dp
203 : IF (do_force) THEN
204 : ALLOCATE (sdabc(ncoa, ncob, ncoc, 3))
205 : sdabc(:, :, :, :) = 0.0_dp
206 : ALLOCATE (sabdc(ncoa, ncob, ncoc, 3))
207 : sabdc(:, :, :, :) = 0.0_dp
208 : CALL overlap3(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
209 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
210 : lc_max(kset), npgfc(kset), zetc(:, kset), rpgfc(:, kset), lc_min(kset), &
211 : rij, dij, rik, dik, rjk, djk, sabc, sdabc, sabdc)
212 : ELSE
213 : CALL overlap3(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
214 : lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
215 : lc_max(kset), npgfc(kset), zetc(:, kset), rpgfc(:, kset), lc_min(kset), &
216 : rij, dij, rik, dik, rjk, djk, sabc)
217 : END IF
218 : ALLOCATE (sabc_contr(nsgfa(iset), nsgfb(jset), nsgfc(kset)))
219 :
220 : CALL abc_contract(sabc_contr, sabc, &
221 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
222 : ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
223 : iabc(sgfa:egfa, sgfb:egfb, sgfc:egfc) = &
224 : sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
225 : IF (do_force) THEN
226 : DO i = 1, 3
227 : CALL abc_contract(sabc_contr, sdabc(:, :, :, i), &
228 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
229 : ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
230 : idabc(sgfa:egfa, sgfb:egfb, sgfc:egfc, i) = &
231 : sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
232 : CALL abc_contract(sabc_contr, sabdc(:, :, :, i), &
233 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
234 : ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
235 : iabdc(sgfa:egfa, sgfb:egfb, sgfc:egfc, i) = &
236 : sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
237 : END DO
238 : END IF
239 :
240 : DEALLOCATE (sabc_contr)
241 : DEALLOCATE (sabc)
242 : END IF
243 : IF (do_force) THEN
244 : DEALLOCATE (sdabc, sabdc)
245 : END IF
246 : END DO
247 : END DO
248 : END DO
249 : IF (do_force) THEN
250 : ! translational invariance
251 : iadbc(:, :, :, :) = -idabc(:, :, :, :) - iabdc(:, :, :, :)
252 : !
253 : ! get the atom indices
254 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, &
255 : iatom=iatom, jatom=jatom, katom=katom)
256 : !
257 : ! contract over i and j to get forces
258 : IF (iatom <= jatom) THEN
259 : irow = iatom
260 : icol = jatom
261 : trans = .FALSE.
262 : ELSE
263 : irow = jatom
264 : icol = iatom
265 : trans = .TRUE.
266 : END IF
267 : IF (iatom == jatom) THEN
268 : fpre = 1.0_dp
269 : ELSE
270 : fpre = 2.0_dp
271 : END IF
272 : ALLOCATE (pmat(ni, nj))
273 : pmat(:, :) = 0.0_dp
274 : DO ispin = 1, nspin
275 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
276 : row=irow, col=icol, BLOCK=pblock, found=found)
277 : IF (found) THEN
278 : IF (trans) THEN
279 : pmat(:, :) = pmat(:, :) + TRANSPOSE(pblock(:, :))
280 : ELSE
281 : pmat(:, :) = pmat(:, :) + pblock(:, :)
282 : END IF
283 : END IF
284 : END DO
285 : DO i = 1, 3
286 : DO j = 1, nk
287 : fi(j, i) = fpre*SUM(pmat(:, :)*idabc(:, :, j, i))
288 : fj(j, i) = fpre*SUM(pmat(:, :)*iadbc(:, :, j, i))
289 : fk(j, i) = fpre*SUM(pmat(:, :)*iabdc(:, :, j, i))
290 : END DO
291 : END DO
292 : DEALLOCATE (pmat)
293 : !
294 : DEALLOCATE (idabc, iadbc, iabdc)
295 : END IF
296 : !
297 : CALL set_o3c_container(o3c_iterator, mepos=mepos, &
298 : integral=iabc, tvec=tvec, force_i=fi, force_j=fj, force_k=fk)
299 :
300 : END DO
301 : !$OMP END PARALLEL
302 0 : CALL o3c_iterator_release(o3c_iterator)
303 :
304 0 : CALL timestop(handle)
305 :
306 0 : END SUBROUTINE calculate_o3c_integrals
307 :
308 : ! **************************************************************************************************
309 : !> \brief Contraction of 3-tensor over indices 1 and 2 (assuming symmetry)
310 : !> t(k) = sum_ij (ijk)*p(ij)
311 : !> \param o3c ...
312 : !> \param matrix_p ...
313 : ! **************************************************************************************************
314 0 : SUBROUTINE contract12_o3c(o3c, matrix_p)
315 : TYPE(o3c_container_type), POINTER :: o3c
316 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
317 :
318 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract12_o3c'
319 :
320 : INTEGER :: handle, iatom, icol, ik, irow, ispin, &
321 : jatom, mepos, nk, nspin, nthread
322 : LOGICAL :: found, ijsymmetric, trans
323 : REAL(KIND=dp) :: fpre
324 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, tvec
325 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: iabc
326 : TYPE(o3c_iterator_type) :: o3c_iterator
327 :
328 0 : CALL timeset(routineN, handle)
329 :
330 0 : nspin = SIZE(matrix_p, 1)
331 0 : CALL get_o3c_container(o3c, ijsymmetric=ijsymmetric)
332 0 : CPASSERT(ijsymmetric)
333 :
334 : nthread = 1
335 0 : !$ nthread = omp_get_max_threads()
336 0 : CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
337 :
338 : !$OMP PARALLEL DEFAULT(NONE) &
339 : !$OMP SHARED (nthread,o3c_iterator,matrix_p,nspin)&
340 0 : !$OMP PRIVATE (mepos,ispin,iatom,jatom,ik,nk,irow,icol,iabc,tvec,found,pblock,trans,fpre)
341 :
342 : mepos = 0
343 : !$ mepos = omp_get_thread_num()
344 :
345 : DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
346 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, iatom=iatom, jatom=jatom, &
347 : integral=iabc, tvec=tvec)
348 : nk = SIZE(tvec, 1)
349 :
350 : IF (iatom <= jatom) THEN
351 : irow = iatom
352 : icol = jatom
353 : trans = .FALSE.
354 : ELSE
355 : irow = jatom
356 : icol = iatom
357 : trans = .TRUE.
358 : END IF
359 : IF (iatom == jatom) THEN
360 : fpre = 1.0_dp
361 : ELSE
362 : fpre = 2.0_dp
363 : END IF
364 :
365 : DO ispin = 1, nspin
366 : CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
367 : row=irow, col=icol, BLOCK=pblock, found=found)
368 : IF (found) THEN
369 : IF (trans) THEN
370 : DO ik = 1, nk
371 : tvec(ik, ispin) = fpre*SUM(TRANSPOSE(pblock(:, :))*iabc(:, :, ik))
372 : END DO
373 : ELSE
374 : DO ik = 1, nk
375 : tvec(ik, ispin) = fpre*SUM(pblock(:, :)*iabc(:, :, ik))
376 : END DO
377 : END IF
378 : END IF
379 : END DO
380 :
381 : END DO
382 : !$OMP END PARALLEL
383 0 : CALL o3c_iterator_release(o3c_iterator)
384 :
385 0 : CALL timestop(handle)
386 :
387 0 : END SUBROUTINE contract12_o3c
388 :
389 : ! **************************************************************************************************
390 : !> \brief Contraction of 3-tensor over index 3
391 : !> h(ij) = h(ij) + sum_k (ijk)*v(k)
392 : !> \param o3c ...
393 : !> \param vec ...
394 : !> \param matrix ...
395 : ! **************************************************************************************************
396 0 : SUBROUTINE contract3_o3c(o3c, vec, matrix)
397 : TYPE(o3c_container_type), POINTER :: o3c
398 : TYPE(o3c_vec_type), DIMENSION(:), POINTER :: vec
399 : TYPE(dbcsr_type) :: matrix
400 :
401 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract3_o3c'
402 :
403 : INTEGER :: handle, iatom, icol, ik, irow, jatom, &
404 : katom, mepos, nk, nthread, s1, s2
405 : LOGICAL :: found, ijsymmetric, trans
406 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
407 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: v
408 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
409 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: iabc
410 : TYPE(o3c_iterator_type) :: o3c_iterator
411 :
412 0 : CALL timeset(routineN, handle)
413 :
414 0 : CALL get_o3c_container(o3c, ijsymmetric=ijsymmetric)
415 0 : CPASSERT(ijsymmetric)
416 :
417 : nthread = 1
418 0 : !$ nthread = omp_get_max_threads()
419 0 : CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
420 :
421 : !$OMP PARALLEL DEFAULT(NONE) &
422 : !$OMP SHARED (nthread,o3c_iterator,vec,matrix)&
423 0 : !$OMP PRIVATE (mepos,iabc,iatom,jatom,katom,irow,icol,trans,pblock,v,found,ik,nk,work,s1,s2)
424 :
425 : mepos = 0
426 : !$ mepos = omp_get_thread_num()
427 :
428 : DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
429 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, iatom=iatom, jatom=jatom, katom=katom, &
430 : integral=iabc)
431 :
432 : CALL get_o3c_vec(vec, katom, v)
433 : nk = SIZE(v)
434 :
435 : IF (iatom <= jatom) THEN
436 : irow = iatom
437 : icol = jatom
438 : trans = .FALSE.
439 : ELSE
440 : irow = jatom
441 : icol = iatom
442 : trans = .TRUE.
443 : END IF
444 :
445 : CALL dbcsr_get_block_p(matrix=matrix, row=irow, col=icol, BLOCK=pblock, found=found)
446 :
447 : IF (found) THEN
448 : s1 = SIZE(pblock, 1); s2 = SIZE(pblock, 2)
449 : ALLOCATE (work(s1, s2))
450 : work(:, :) = 0.0_dp
451 :
452 : IF (trans) THEN
453 : DO ik = 1, nk
454 : CALL daxpy(s1*s2, v(ik), TRANSPOSE(iabc(:, :, ik)), 1, work(:, :), 1)
455 : END DO
456 : ELSE
457 : DO ik = 1, nk
458 : CALL daxpy(s1*s2, v(ik), iabc(:, :, ik), 1, work(:, :), 1)
459 : END DO
460 : END IF
461 :
462 : ! Multiple threads with same irow, icol but different katom (same even in PBCs) can try
463 : ! to access the dbcsr block at the same time. Prevent that by CRITICAL section but keep
464 : ! computations before hand in order to retain speed
465 :
466 : !$OMP CRITICAL
467 : CALL dbcsr_get_block_p(matrix=matrix, row=irow, col=icol, BLOCK=pblock, found=found)
468 : CALL daxpy(s1*s2, 1.0_dp, work(:, :), 1, pblock(:, :), 1)
469 : !$OMP END CRITICAL
470 :
471 : DEALLOCATE (work)
472 : END IF
473 :
474 : END DO
475 : !$OMP END PARALLEL
476 0 : CALL o3c_iterator_release(o3c_iterator)
477 :
478 0 : CALL timestop(handle)
479 :
480 0 : END SUBROUTINE contract3_o3c
481 :
482 : END MODULE qs_o3c_methods
|
