Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Utility method to build 3-center integrals for small cell GW
10 : ! **************************************************************************************************
11 : MODULE gw_integrals
12 : USE OMP_LIB, ONLY: omp_get_thread_num
13 : USE ai_contraction_sphi, ONLY: abc_contract_xsmm
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind_set
16 : USE basis_set_types, ONLY: get_gto_basis_set,&
17 : gto_basis_set_p_type,&
18 : gto_basis_set_type
19 : USE cell_types, ONLY: cell_type,&
20 : get_cell,&
21 : pbc
22 : USE cp_array_utils, ONLY: cp_2d_r_p_type
23 : USE cp_files, ONLY: close_file,&
24 : open_file
25 : USE gamma, ONLY: init_md_ftable
26 : USE input_constants, ONLY: do_potential_coulomb,&
27 : do_potential_id,&
28 : do_potential_short,&
29 : do_potential_truncated
30 : USE kinds, ONLY: dp
31 : USE libint_2c_3c, ONLY: cutoff_screen_factor,&
32 : eri_3center,&
33 : libint_potential_type
34 : USE libint_wrapper, ONLY: cp_libint_cleanup_3eri,&
35 : cp_libint_init_3eri,&
36 : cp_libint_set_contrdepth,&
37 : cp_libint_t
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE orbital_pointers, ONLY: ncoset
40 : USE particle_types, ONLY: particle_type
41 : USE qs_environment_types, ONLY: get_qs_env,&
42 : qs_environment_type
43 : USE t_c_g0, ONLY: get_lmax_init,&
44 : init
45 :
46 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 :
51 : PRIVATE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_integrals'
54 :
55 : PUBLIC :: build_3c_integral_block, build_3c_integral_block_ctx, &
56 : gw_3c_ctx_type, gw_3c_ctx_create, gw_3c_ctx_release, &
57 : gw_3c_ws_type, gw_3c_ws_create, gw_3c_ws_release
58 :
59 : ! **************************************************************************************************
60 : !> \brief Shared read-only context for repeated 3-center integral block builds: screening
61 : !> parameters, basis maxima, contracted sphi tables, and the one-time gamma /
62 : !> truncated-Coulomb table initializations. Create and release OUTSIDE any OMP parallel
63 : !> region; creation is MPI-collective when the potential is truncated.
64 : ! **************************************************************************************************
65 : TYPE gw_3c_ctx_type
66 : TYPE(libint_potential_type) :: potential_parameter = libint_potential_type()
67 : INTEGER :: op_ij = do_potential_id, &
68 : op_jk = do_potential_id
69 : REAL(KIND=dp) :: dr_ij = 0.0_dp, dr_jk = 0.0_dp, &
70 : dr_ik = 0.0_dp
71 : INTEGER :: maxli = 0, maxlj = 0, maxlk = 0, &
72 : max_am = 0, m_max = 0
73 : INTEGER :: max_ncoi = 0, max_ncoj = 0, max_ncok = 0
74 : INTEGER :: max_nsgfi = 0, max_nsgfj = 0, &
75 : max_nsgfk = 0, max_nset = 0, natom = 0
76 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi => NULL(), tspj => NULL(), &
77 : spk => NULL()
78 : TYPE(gto_basis_set_p_type), DIMENSION(:), ALLOCATABLE :: basis_i, basis_j, basis_k
79 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
80 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set => NULL()
81 : TYPE(cell_type), POINTER :: cell => NULL()
82 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat = 0.0_dp
83 : END TYPE gw_3c_ctx_type
84 :
85 : ! **************************************************************************************************
86 : !> \brief Per-thread workspace for 3-center integral block builds: libint object + contraction
87 : !> buffers. Each thread creates its own (inside the parallel region is fine).
88 : ! **************************************************************************************************
89 : TYPE gw_3c_ws_type
90 : TYPE(cp_libint_t), ALLOCATABLE :: lib
91 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpp_buffer, ccp_buffer
92 : END TYPE gw_3c_ws_type
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief Builds the shared 3c-integral context: screening radii, basis maxima, contracted sphi
98 : !> tables, and the one-time gamma / truncated-Coulomb table initializations.
99 : !> \param ctx ...
100 : !> \param qs_env ...
101 : !> \param potential_parameter ...
102 : !> \param basis_j ...
103 : !> \param basis_k ...
104 : !> \param basis_i ...
105 : ! **************************************************************************************************
106 181664 : SUBROUTINE gw_3c_ctx_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i)
107 :
108 : TYPE(gw_3c_ctx_type), INTENT(OUT) :: ctx
109 : TYPE(qs_environment_type), POINTER :: qs_env
110 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
111 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_j, basis_k, basis_i
112 :
113 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_3c_ctx_create'
114 :
115 : INTEGER :: egfi, handle, ibasis, ilist, imax, iset, &
116 : jset, kset, l, nbasis, ncoi, npgf_l, &
117 : sgfi, unit_id
118 12976 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmin_i, npgfi, npgfj, npgfk, &
119 12976 : nsgfi, nsgfj, nsgfk
120 12976 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
121 : TYPE(gto_basis_set_type), POINTER :: basis_set
122 : TYPE(mp_para_env_type), POINTER :: para_env
123 :
124 12976 : CALL timeset(routineN, handle)
125 :
126 12976 : ctx%potential_parameter = potential_parameter
127 12976 : ctx%op_ij = potential_parameter%potential_type
128 12976 : ctx%op_jk = do_potential_id
129 :
130 12976 : IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_ij == do_potential_short) THEN
131 12976 : ctx%dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
132 12976 : ctx%dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
133 0 : ELSE IF (ctx%op_ij == do_potential_coulomb) THEN
134 0 : ctx%dr_ij = 1000000.0_dp
135 0 : ctx%dr_ik = 1000000.0_dp
136 : END IF
137 :
138 12976 : NULLIFY (atomic_kind_set, para_env)
139 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=ctx%cell, natom=ctx%natom, &
140 12976 : para_env=para_env, particle_set=ctx%particle_set)
141 12976 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=ctx%kind_of)
142 12976 : CALL get_cell(cell=ctx%cell, h=ctx%hmat)
143 :
144 51778 : ctx%basis_i = basis_i
145 51778 : ctx%basis_j = basis_j
146 51778 : ctx%basis_k = basis_k
147 :
148 : ! max l per basis for libint; max nset/nco/nsgf for the LIBXSMM contraction buffers
149 12976 : nbasis = SIZE(basis_i)
150 38802 : DO ibasis = 1, nbasis
151 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
152 25826 : lmax=lmax_i, lmin=lmin_i, nset=iset, nsgf_set=nsgfi, npgf=npgfi)
153 25826 : ctx%maxli = MAX(ctx%maxli, imax)
154 25826 : ctx%max_nset = MAX(ctx%max_nset, iset)
155 65732 : ctx%max_nsgfi = MAX(ctx%max_nsgfi, MAXVAL(nsgfi))
156 65732 : ctx%max_ncoi = MAX(ctx%max_ncoi, MAXVAL(npgfi)*ncoset(ctx%maxli))
157 103620 : DO l = 0, imax
158 38992 : npgf_l = 0
159 107900 : DO jset = 1, iset
160 107900 : IF (lmin_i(jset) == l .AND. lmax_i(jset) == l) npgf_l = npgf_l + npgfi(jset)
161 : END DO
162 64818 : ctx%max_ncoi = MAX(ctx%max_ncoi, npgf_l*ncoset(l))
163 : END DO
164 : END DO
165 38802 : DO ibasis = 1, nbasis
166 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
167 25826 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
168 25826 : ctx%maxlj = MAX(ctx%maxlj, imax)
169 25826 : ctx%max_nset = MAX(ctx%max_nset, jset)
170 77610 : ctx%max_nsgfj = MAX(ctx%max_nsgfj, MAXVAL(nsgfj))
171 116412 : ctx%max_ncoj = MAX(ctx%max_ncoj, MAXVAL(npgfj)*ncoset(ctx%maxlj))
172 : END DO
173 38802 : DO ibasis = 1, nbasis
174 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
175 25826 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
176 25826 : ctx%maxlk = MAX(ctx%maxlk, imax)
177 25826 : ctx%max_nset = MAX(ctx%max_nset, kset)
178 77610 : ctx%max_nsgfk = MAX(ctx%max_nsgfk, MAXVAL(nsgfk))
179 116412 : ctx%max_ncok = MAX(ctx%max_ncok, MAXVAL(npgfk)*ncoset(ctx%maxlk))
180 : END DO
181 12976 : ctx%m_max = ctx%maxli + ctx%maxlj + ctx%maxlk
182 12976 : ctx%max_am = MAX(ctx%maxli, ctx%maxlj, ctx%maxlk)
183 :
184 : ! contiguous (and for j transposed) sphi copies, shared read-only across threads
185 : ALLOCATE (ctx%spi(ctx%max_nset, nbasis), ctx%tspj(ctx%max_nset, nbasis), &
186 442720 : ctx%spk(ctx%max_nset, nbasis))
187 38802 : DO ibasis = 1, nbasis
188 117296 : DO iset = 1, ctx%max_nset
189 78494 : NULLIFY (ctx%spi(iset, ibasis)%array)
190 78494 : NULLIFY (ctx%tspj(iset, ibasis)%array)
191 104320 : NULLIFY (ctx%spk(iset, ibasis)%array)
192 : END DO
193 : END DO
194 51904 : DO ilist = 1, 3
195 129382 : DO ibasis = 1, nbasis
196 77478 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
197 77478 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
198 77478 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
199 259880 : DO iset = 1, basis_set%nset
200 143474 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
201 143474 : sgfi = basis_set%first_sgf(1, iset)
202 143474 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
203 220952 : IF (ilist == 1) THEN
204 159624 : ALLOCATE (ctx%spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
205 324338 : ctx%spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
206 103568 : ELSE IF (ilist == 2) THEN
207 207136 : ALLOCATE (ctx%tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
208 1646452 : ctx%tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
209 : ELSE
210 207136 : ALLOCATE (ctx%spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
211 1402196 : ctx%spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
212 : END IF
213 : END DO
214 : END DO
215 : END DO
216 :
217 : ! one-time table inits; the truncated-Coulomb init reads a file + bcasts => MPI-collective,
218 : ! must happen here and never inside the per-block path or an OMP region
219 12976 : IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_jk == do_potential_truncated) THEN
220 12976 : IF (ctx%m_max > get_lmax_init()) THEN
221 0 : IF (para_env%mepos == 0) THEN
222 0 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
223 : END IF
224 0 : CALL init(ctx%m_max, unit_id, para_env%mepos, para_env)
225 0 : IF (para_env%mepos == 0) THEN
226 0 : CALL close_file(unit_id)
227 : END IF
228 : END IF
229 : END IF
230 12976 : CALL init_md_ftable(nmax=ctx%m_max)
231 :
232 12976 : CALL timestop(handle)
233 :
234 25952 : END SUBROUTINE gw_3c_ctx_create
235 :
236 : ! **************************************************************************************************
237 : !> \brief Releases the shared 3c-integral context.
238 : !> \param ctx ...
239 : ! **************************************************************************************************
240 12976 : SUBROUTINE gw_3c_ctx_release(ctx)
241 :
242 : TYPE(gw_3c_ctx_type), INTENT(INOUT) :: ctx
243 :
244 : INTEGER :: ibasis, iset
245 :
246 52726 : DO iset = 1, SIZE(ctx%spi, 1)
247 131220 : DO ibasis = 1, SIZE(ctx%spi, 2)
248 78494 : IF (ASSOCIATED(ctx%spi(iset, ibasis)%array)) DEALLOCATE (ctx%spi(iset, ibasis)%array)
249 78494 : IF (ASSOCIATED(ctx%tspj(iset, ibasis)%array)) DEALLOCATE (ctx%tspj(iset, ibasis)%array)
250 118244 : IF (ASSOCIATED(ctx%spk(iset, ibasis)%array)) DEALLOCATE (ctx%spk(iset, ibasis)%array)
251 : END DO
252 : END DO
253 12976 : DEALLOCATE (ctx%spi, ctx%tspj, ctx%spk)
254 12976 : NULLIFY (ctx%spi, ctx%tspj, ctx%spk, ctx%particle_set, ctx%cell)
255 12976 : IF (ALLOCATED(ctx%kind_of)) DEALLOCATE (ctx%kind_of)
256 12976 : IF (ALLOCATED(ctx%basis_i)) DEALLOCATE (ctx%basis_i)
257 12976 : IF (ALLOCATED(ctx%basis_j)) DEALLOCATE (ctx%basis_j)
258 12976 : IF (ALLOCATED(ctx%basis_k)) DEALLOCATE (ctx%basis_k)
259 :
260 12976 : END SUBROUTINE gw_3c_ctx_release
261 :
262 : ! **************************************************************************************************
263 : !> \brief Creates a per-thread 3c workspace: libint object + LIBXSMM contraction buffers.
264 : !> \param ws ...
265 : !> \param ctx ...
266 : ! **************************************************************************************************
267 12985 : SUBROUTINE gw_3c_ws_create(ws, ctx)
268 :
269 : TYPE(gw_3c_ws_type), INTENT(OUT) :: ws
270 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
271 :
272 12985 : ALLOCATE (ws%lib)
273 12985 : CALL cp_libint_init_3eri(ws%lib, ctx%max_am)
274 12985 : CALL cp_libint_set_contrdepth(ws%lib, 1)
275 0 : ALLOCATE (ws%cpp_buffer(ctx%max_nsgfj*ctx%max_ncok), &
276 64925 : ws%ccp_buffer(ctx%max_nsgfj*ctx%max_nsgfk*ctx%max_ncoi))
277 :
278 12985 : END SUBROUTINE gw_3c_ws_create
279 :
280 : ! **************************************************************************************************
281 : !> \brief Releases a per-thread 3c workspace.
282 : !> \param ws ...
283 : ! **************************************************************************************************
284 12985 : SUBROUTINE gw_3c_ws_release(ws)
285 :
286 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
287 :
288 12985 : CALL cp_libint_cleanup_3eri(ws%lib)
289 12985 : DEALLOCATE (ws%lib)
290 12985 : DEALLOCATE (ws%cpp_buffer, ws%ccp_buffer)
291 :
292 12985 : END SUBROUTINE gw_3c_ws_release
293 :
294 : ! **************************************************************************************************
295 : !> \brief Computes the 3c integral block (mu(atom_j) nu(atom_k) | P(atom_i)) for ONE atom triple,
296 : !> accumulating into the caller-zeroed int_3c at the given block offsets.
297 : !> Thread-safe: reads the frozen ctx + read-only module tables, mutates only its arguments
298 : !> and the per-thread ws.
299 : !> \param int_3c pre-zeroed target; the triple's contribution is accumulated in place
300 : !> \param ctx shared context from gw_3c_ctx_create
301 : !> \param ws per-thread workspace from gw_3c_ws_create
302 : !> \param atom_j ...
303 : !> \param atom_k ...
304 : !> \param atom_i ...
305 : !> \param cell_j ...
306 : !> \param cell_k ...
307 : !> \param cell_i ...
308 : !> \param j_offset block offset of atom_j's first sgf in int_3c dim 1 (default 0)
309 : !> \param k_offset block offset of atom_k's first sgf in int_3c dim 2 (default 0)
310 : !> \param i_offset block offset of atom_i's first RI sgf in int_3c dim 3 (default 0)
311 : !> \param screened .TRUE. if the kind-radius screens killed the whole triple (int_3c untouched)
312 : ! **************************************************************************************************
313 13592 : SUBROUTINE build_3c_integral_block_ctx(int_3c, ctx, ws, atom_j, atom_k, atom_i, &
314 : cell_j, cell_k, cell_i, &
315 : j_offset, k_offset, i_offset, screened)
316 :
317 : REAL(KIND=dp), DIMENSION(:, :, :) :: int_3c
318 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
319 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
320 : INTEGER, INTENT(IN) :: atom_j, atom_k, atom_i
321 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_j, cell_k, cell_i
322 : INTEGER, INTENT(IN), OPTIONAL :: j_offset, k_offset, i_offset
323 : LOGICAL, INTENT(OUT), OPTIONAL :: screened
324 :
325 : INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
326 : block_start_k, ikind, iset, jkind, jset, kkind, kset, my_i_offset, my_j_offset, &
327 : my_k_offset, ncoi, ncoj, ncok, nseti, nsetj, nsetk, sgfi, sgfj, sgfk
328 : INTEGER, DIMENSION(3) :: my_cell_i, my_cell_j, my_cell_k
329 13592 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
330 13592 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
331 13592 : nsgfj, nsgfk
332 13592 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
333 : REAL(KIND=dp) :: dij, dik, djk, kind_radius_i, &
334 : kind_radius_j, kind_radius_k, sijk_ext
335 13592 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sijk, sijk_contr
336 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
337 13592 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
338 13592 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, zeti, zetj, zetk
339 :
340 570 : IF (PRESENT(screened)) screened = .FALSE.
341 :
342 13592 : my_cell_i(1:3) = 0
343 13592 : IF (PRESENT(cell_i)) my_cell_i(1:3) = cell_i(1:3)
344 13592 : my_cell_j(1:3) = 0
345 13592 : IF (PRESENT(cell_j)) my_cell_j(1:3) = cell_j(1:3)
346 13592 : my_cell_k(1:3) = 0
347 13592 : IF (PRESENT(cell_k)) my_cell_k(1:3) = cell_k(1:3)
348 13592 : my_i_offset = 0
349 13592 : IF (PRESENT(i_offset)) my_i_offset = i_offset
350 13592 : my_j_offset = 0
351 13592 : IF (PRESENT(j_offset)) my_j_offset = j_offset
352 13592 : my_k_offset = 0
353 13592 : IF (PRESENT(k_offset)) my_k_offset = k_offset
354 :
355 271840 : ri = pbc(ctx%particle_set(atom_i)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_i, dp))
356 271840 : rj = pbc(ctx%particle_set(atom_j)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_j, dp))
357 271840 : rk = pbc(ctx%particle_set(atom_k)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_k, dp))
358 :
359 54368 : rjk(1:3) = rk(1:3) - rj(1:3)
360 54368 : rij(1:3) = rj(1:3) - ri(1:3)
361 54368 : rik(1:3) = rk(1:3) - ri(1:3)
362 :
363 54368 : djk = NORM2(rjk)
364 54368 : dij = NORM2(rij)
365 54368 : dik = NORM2(rik)
366 :
367 13592 : ikind = ctx%kind_of(atom_i)
368 13592 : jkind = ctx%kind_of(atom_j)
369 13592 : kkind = ctx%kind_of(atom_k)
370 :
371 : CALL get_gto_basis_set(ctx%basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, &
372 : lmax=lmax_i, lmin=lmin_i, npgf=npgfi, nset=nseti, &
373 : nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
374 13592 : zet=zeti, kind_radius=kind_radius_i)
375 : CALL get_gto_basis_set(ctx%basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, &
376 : lmax=lmax_j, lmin=lmin_j, npgf=npgfj, nset=nsetj, &
377 : nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
378 13592 : zet=zetj, kind_radius=kind_radius_j)
379 : CALL get_gto_basis_set(ctx%basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, &
380 : lmax=lmax_k, lmin=lmin_k, npgf=npgfk, nset=nsetk, &
381 : nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
382 13592 : zet=zetk, kind_radius=kind_radius_k)
383 :
384 : IF (kind_radius_j + kind_radius_i + ctx%dr_ij < dij .OR. &
385 13592 : kind_radius_j + kind_radius_k + ctx%dr_jk < djk .OR. &
386 : kind_radius_k + kind_radius_i + ctx%dr_ik < dik) THEN
387 0 : IF (PRESENT(screened)) screened = .TRUE.
388 0 : RETURN
389 : END IF
390 :
391 36156 : DO iset = 1, nseti
392 85507 : DO jset = 1, nsetj
393 49351 : IF (set_radius_j(jset) + set_radius_i(iset) + ctx%dr_ij < dij) CYCLE
394 177246 : DO kset = 1, nsetk
395 109075 : IF (set_radius_j(jset) + set_radius_k(kset) + ctx%dr_jk < djk) CYCLE
396 100975 : IF (set_radius_k(kset) + set_radius_i(iset) + ctx%dr_ik < dik) CYCLE
397 :
398 95433 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
399 95433 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
400 95433 : ncok = npgfk(kset)*ncoset(lmax_k(kset))
401 :
402 95433 : sgfi = first_sgf_i(1, iset)
403 95433 : sgfj = first_sgf_j(1, jset)
404 95433 : sgfk = first_sgf_k(1, kset)
405 :
406 95433 : IF (ncoj*ncok*ncoi <= 0) CYCLE
407 477165 : ALLOCATE (sijk(ncoj, ncok, ncoi))
408 95433 : sijk(:, :, :) = 0.0_dp
409 :
410 : CALL eri_3center(sijk, &
411 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
412 : rpgf_j(:, jset), rj, &
413 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), &
414 : rpgf_k(:, kset), rk, &
415 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
416 : rpgf_i(:, iset), ri, &
417 : djk, dij, dik, ws%lib, ctx%potential_parameter, &
418 95433 : int_abc_ext=sijk_ext)
419 :
420 477165 : ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
421 : CALL abc_contract_xsmm(sijk_contr, sijk, ctx%tspj(jset, jkind)%array, &
422 : ctx%spk(kset, kkind)%array, ctx%spi(iset, ikind)%array, &
423 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
424 95433 : nsgfi(iset), ws%cpp_buffer, ws%ccp_buffer)
425 95433 : DEALLOCATE (sijk)
426 :
427 95433 : block_start_j = sgfj + my_j_offset
428 95433 : block_end_j = sgfj + nsgfj(jset) - 1 + my_j_offset
429 95433 : block_start_k = sgfk + my_k_offset
430 95433 : block_end_k = sgfk + nsgfk(kset) - 1 + my_k_offset
431 95433 : block_start_i = sgfi + my_i_offset
432 95433 : block_end_i = sgfi + nsgfi(iset) - 1 + my_i_offset
433 :
434 : int_3c(block_start_j:block_end_j, &
435 : block_start_k:block_end_k, &
436 : block_start_i:block_end_i) = &
437 : int_3c(block_start_j:block_end_j, &
438 : block_start_k:block_end_k, &
439 : block_start_i:block_end_i) + &
440 1239115 : sijk_contr(:, :, :)
441 158426 : DEALLOCATE (sijk_contr)
442 :
443 : END DO
444 : END DO
445 : END DO
446 :
447 40776 : END SUBROUTINE build_3c_integral_block_ctx
448 :
449 : ! **************************************************************************************************
450 : !> \brief ...
451 : !> \param int_3c ...
452 : !> \param qs_env ...
453 : !> \param potential_parameter ...
454 : !> \param basis_j ...
455 : !> \param basis_k ...
456 : !> \param basis_i ...
457 : !> \param cell_j ...
458 : !> \param cell_k ...
459 : !> \param cell_i ...
460 : !> \param atom_j ...
461 : !> \param atom_k ...
462 : !> \param atom_i ...
463 : !> \param j_bf_start_from_atom ...
464 : !> \param k_bf_start_from_atom ...
465 : !> \param i_bf_start_from_atom ...
466 : ! **************************************************************************************************
467 38766 : SUBROUTINE build_3c_integral_block(int_3c, qs_env, potential_parameter, &
468 12922 : basis_j, basis_k, basis_i, &
469 : cell_j, cell_k, cell_i, atom_j, atom_k, atom_i, &
470 12922 : j_bf_start_from_atom, k_bf_start_from_atom, &
471 12922 : i_bf_start_from_atom)
472 :
473 : REAL(KIND=dp), DIMENSION(:, :, :) :: int_3c
474 : TYPE(qs_environment_type), POINTER :: qs_env
475 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
476 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_j, basis_k, basis_i
477 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_j, cell_k, cell_i
478 : INTEGER, INTENT(IN), OPTIONAL :: atom_j, atom_k, atom_i
479 : INTEGER, DIMENSION(:), OPTIONAL :: j_bf_start_from_atom, &
480 : k_bf_start_from_atom, &
481 : i_bf_start_from_atom
482 :
483 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integral_block'
484 :
485 : INTEGER :: at_i, at_j, at_k, handle, my_i_offset, &
486 : my_j_offset, my_k_offset
487 167986 : TYPE(gw_3c_ctx_type) :: ctx
488 12922 : TYPE(gw_3c_ws_type) :: ws
489 :
490 12922 : CALL timeset(routineN, handle)
491 :
492 12922 : CALL gw_3c_ctx_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i)
493 12922 : CALL gw_3c_ws_create(ws, ctx)
494 :
495 823215 : int_3c(:, :, :) = 0.0_dp
496 :
497 : ! loop over all RI atoms
498 51303 : DO at_i = 1, ctx%natom
499 : ! loop over all AO atoms
500 165676 : DO at_j = 1, ctx%natom
501 : ! loop over all AO atoms
502 494333 : DO at_k = 1, ctx%natom
503 :
504 341579 : IF (PRESENT(atom_i)) THEN
505 340779 : IF (at_i /= atom_i) CYCLE
506 : END IF
507 114773 : IF (PRESENT(atom_j)) THEN
508 114773 : IF (at_j /= atom_j) CYCLE
509 : END IF
510 38581 : IF (PRESENT(atom_k)) THEN
511 38581 : IF (at_k /= atom_k) CYCLE
512 : END IF
513 :
514 13022 : IF (PRESENT(atom_j)) THEN
515 13022 : my_j_offset = 0
516 : ELSE
517 0 : CPASSERT(PRESENT(j_bf_start_from_atom))
518 0 : my_j_offset = j_bf_start_from_atom(at_j) - 1
519 : END IF
520 13022 : IF (PRESENT(atom_k)) THEN
521 13022 : my_k_offset = 0
522 : ELSE
523 0 : CPASSERT(PRESENT(k_bf_start_from_atom))
524 0 : my_k_offset = k_bf_start_from_atom(at_k) - 1
525 : END IF
526 13022 : IF (PRESENT(atom_i)) THEN
527 12822 : my_i_offset = 0
528 : ELSE
529 200 : CPASSERT(PRESENT(i_bf_start_from_atom))
530 200 : my_i_offset = i_bf_start_from_atom(at_i) - 1
531 : END IF
532 :
533 : CALL build_3c_integral_block_ctx(int_3c, ctx, ws, at_j, at_k, at_i, &
534 : cell_j=cell_j, cell_k=cell_k, cell_i=cell_i, &
535 : j_offset=my_j_offset, k_offset=my_k_offset, &
536 455952 : i_offset=my_i_offset)
537 :
538 : END DO ! atom_k (AO)
539 : END DO ! atom_j (AO)
540 : END DO ! atom_i (RI)
541 :
542 12922 : CALL gw_3c_ws_release(ws)
543 12922 : CALL gw_3c_ctx_release(ctx)
544 :
545 12922 : CALL timestop(handle)
546 :
547 12922 : END SUBROUTINE build_3c_integral_block
548 :
549 0 : END MODULE gw_integrals
|
