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) :: 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 180040 : 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, nbasis, ncoi, sgfi, unit_id
117 12860 : INTEGER, DIMENSION(:), POINTER :: npgfi, npgfj, npgfk, nsgfi, nsgfj, nsgfk
118 12860 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
119 : TYPE(gto_basis_set_type), POINTER :: basis_set
120 : TYPE(mp_para_env_type), POINTER :: para_env
121 :
122 12860 : CALL timeset(routineN, handle)
123 :
124 12860 : ctx%potential_parameter = potential_parameter
125 12860 : ctx%op_ij = potential_parameter%potential_type
126 12860 : ctx%op_jk = do_potential_id
127 :
128 12860 : IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_ij == do_potential_short) THEN
129 12860 : ctx%dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
130 12860 : ctx%dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
131 0 : ELSE IF (ctx%op_ij == do_potential_coulomb) THEN
132 0 : ctx%dr_ij = 1000000.0_dp
133 0 : ctx%dr_ik = 1000000.0_dp
134 : END IF
135 :
136 12860 : NULLIFY (atomic_kind_set, para_env)
137 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=ctx%cell, natom=ctx%natom, &
138 12860 : para_env=para_env, particle_set=ctx%particle_set)
139 12860 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=ctx%kind_of)
140 12860 : CALL get_cell(cell=ctx%cell, h=ctx%hmat)
141 :
142 51410 : ctx%basis_i = basis_i
143 51410 : ctx%basis_j = basis_j
144 51410 : ctx%basis_k = basis_k
145 :
146 : ! max l per basis for libint; max nset/nco/nsgf for the LIBXSMM contraction buffers
147 12860 : nbasis = SIZE(basis_i)
148 38550 : DO ibasis = 1, nbasis
149 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
150 25690 : nset=iset, nsgf_set=nsgfi, npgf=npgfi)
151 25690 : ctx%maxli = MAX(ctx%maxli, imax)
152 25690 : ctx%max_nset = MAX(ctx%max_nset, iset)
153 64370 : ctx%max_nsgfi = MAX(ctx%max_nsgfi, MAXVAL(nsgfi))
154 102920 : ctx%max_ncoi = MAX(ctx%max_ncoi, MAXVAL(npgfi)*ncoset(ctx%maxli))
155 : END DO
156 38550 : DO ibasis = 1, nbasis
157 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
158 25690 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
159 25690 : ctx%maxlj = MAX(ctx%maxlj, imax)
160 25690 : ctx%max_nset = MAX(ctx%max_nset, jset)
161 77034 : ctx%max_nsgfj = MAX(ctx%max_nsgfj, MAXVAL(nsgfj))
162 115584 : ctx%max_ncoj = MAX(ctx%max_ncoj, MAXVAL(npgfj)*ncoset(ctx%maxlj))
163 : END DO
164 38550 : DO ibasis = 1, nbasis
165 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
166 25690 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
167 25690 : ctx%maxlk = MAX(ctx%maxlk, imax)
168 25690 : ctx%max_nset = MAX(ctx%max_nset, kset)
169 77034 : ctx%max_nsgfk = MAX(ctx%max_nsgfk, MAXVAL(nsgfk))
170 115584 : ctx%max_ncok = MAX(ctx%max_ncok, MAXVAL(npgfk)*ncoset(ctx%maxlk))
171 : END DO
172 12860 : ctx%m_max = ctx%maxli + ctx%maxlj + ctx%maxlk
173 12860 : ctx%max_am = MAX(ctx%maxli, ctx%maxlj, ctx%maxlk)
174 :
175 : ! contiguous (and for j transposed) sphi copies, shared read-only across threads
176 : ALLOCATE (ctx%spi(ctx%max_nset, nbasis), ctx%tspj(ctx%max_nset, nbasis), &
177 437156 : ctx%spk(ctx%max_nset, nbasis))
178 38550 : DO ibasis = 1, nbasis
179 115712 : DO iset = 1, ctx%max_nset
180 77162 : NULLIFY (ctx%spi(iset, ibasis)%array)
181 77162 : NULLIFY (ctx%tspj(iset, ibasis)%array)
182 102852 : NULLIFY (ctx%spk(iset, ibasis)%array)
183 : END DO
184 : END DO
185 51440 : DO ilist = 1, 3
186 128510 : DO ibasis = 1, nbasis
187 77070 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
188 77070 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
189 77070 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
190 257018 : DO iset = 1, basis_set%nset
191 141368 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
192 141368 : sgfi = basis_set%first_sgf(1, iset)
193 141368 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
194 218438 : IF (ilist == 1) THEN
195 154720 : ALLOCATE (ctx%spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
196 284582 : ctx%spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
197 102688 : ELSE IF (ilist == 2) THEN
198 205376 : ALLOCATE (ctx%tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
199 1642372 : ctx%tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
200 : ELSE
201 205376 : ALLOCATE (ctx%spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
202 1398470 : ctx%spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
203 : END IF
204 : END DO
205 : END DO
206 : END DO
207 :
208 : ! one-time table inits; the truncated-Coulomb init reads a file + bcasts => MPI-collective,
209 : ! must happen here and never inside the per-block path or an OMP region
210 12860 : IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_jk == do_potential_truncated) THEN
211 12860 : IF (ctx%m_max > get_lmax_init()) THEN
212 16 : IF (para_env%mepos == 0) THEN
213 8 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
214 : END IF
215 16 : CALL init(ctx%m_max, unit_id, para_env%mepos, para_env)
216 16 : IF (para_env%mepos == 0) THEN
217 8 : CALL close_file(unit_id)
218 : END IF
219 : END IF
220 : END IF
221 12860 : CALL init_md_ftable(nmax=ctx%m_max)
222 :
223 12860 : CALL timestop(handle)
224 :
225 25720 : END SUBROUTINE gw_3c_ctx_create
226 :
227 : ! **************************************************************************************************
228 : !> \brief Releases the shared 3c-integral context.
229 : !> \param ctx ...
230 : ! **************************************************************************************************
231 12860 : SUBROUTINE gw_3c_ctx_release(ctx)
232 :
233 : TYPE(gw_3c_ctx_type), INTENT(INOUT) :: ctx
234 :
235 : INTEGER :: ibasis, iset
236 :
237 51516 : DO iset = 1, SIZE(ctx%spi, 1)
238 128678 : DO ibasis = 1, SIZE(ctx%spi, 2)
239 77162 : IF (ASSOCIATED(ctx%spi(iset, ibasis)%array)) DEALLOCATE (ctx%spi(iset, ibasis)%array)
240 77162 : IF (ASSOCIATED(ctx%tspj(iset, ibasis)%array)) DEALLOCATE (ctx%tspj(iset, ibasis)%array)
241 115818 : IF (ASSOCIATED(ctx%spk(iset, ibasis)%array)) DEALLOCATE (ctx%spk(iset, ibasis)%array)
242 : END DO
243 : END DO
244 12860 : DEALLOCATE (ctx%spi, ctx%tspj, ctx%spk)
245 12860 : NULLIFY (ctx%spi, ctx%tspj, ctx%spk, ctx%particle_set, ctx%cell)
246 12860 : IF (ALLOCATED(ctx%kind_of)) DEALLOCATE (ctx%kind_of)
247 12860 : IF (ALLOCATED(ctx%basis_i)) DEALLOCATE (ctx%basis_i)
248 12860 : IF (ALLOCATED(ctx%basis_j)) DEALLOCATE (ctx%basis_j)
249 12860 : IF (ALLOCATED(ctx%basis_k)) DEALLOCATE (ctx%basis_k)
250 :
251 12860 : END SUBROUTINE gw_3c_ctx_release
252 :
253 : ! **************************************************************************************************
254 : !> \brief Creates a per-thread 3c workspace: libint object + LIBXSMM contraction buffers.
255 : !> \param ws ...
256 : !> \param ctx ...
257 : ! **************************************************************************************************
258 12864 : SUBROUTINE gw_3c_ws_create(ws, ctx)
259 :
260 : TYPE(gw_3c_ws_type), INTENT(OUT) :: ws
261 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
262 :
263 12864 : CALL cp_libint_init_3eri(ws%lib, ctx%max_am)
264 12864 : CALL cp_libint_set_contrdepth(ws%lib, 1)
265 0 : ALLOCATE (ws%cpp_buffer(ctx%max_nsgfj*ctx%max_ncok), &
266 64320 : ws%ccp_buffer(ctx%max_nsgfj*ctx%max_nsgfk*ctx%max_ncoi))
267 :
268 12864 : END SUBROUTINE gw_3c_ws_create
269 :
270 : ! **************************************************************************************************
271 : !> \brief Releases a per-thread 3c workspace.
272 : !> \param ws ...
273 : ! **************************************************************************************************
274 12864 : SUBROUTINE gw_3c_ws_release(ws)
275 :
276 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
277 :
278 12864 : CALL cp_libint_cleanup_3eri(ws%lib)
279 12864 : DEALLOCATE (ws%cpp_buffer, ws%ccp_buffer)
280 :
281 12864 : END SUBROUTINE gw_3c_ws_release
282 :
283 : ! **************************************************************************************************
284 : !> \brief Computes the 3c integral block (mu(atom_j) nu(atom_k) | P(atom_i)) for ONE atom triple,
285 : !> accumulating into the caller-zeroed int_3c at the given block offsets.
286 : !> Thread-safe: reads the frozen ctx + read-only module tables, mutates only its arguments
287 : !> and the per-thread ws.
288 : !> \param int_3c pre-zeroed target; the triple's contribution is accumulated in place
289 : !> \param ctx shared context from gw_3c_ctx_create
290 : !> \param ws per-thread workspace from gw_3c_ws_create
291 : !> \param atom_j ...
292 : !> \param atom_k ...
293 : !> \param atom_i ...
294 : !> \param cell_j ...
295 : !> \param cell_k ...
296 : !> \param cell_i ...
297 : !> \param j_offset block offset of atom_j's first sgf in int_3c dim 1 (default 0)
298 : !> \param k_offset block offset of atom_k's first sgf in int_3c dim 2 (default 0)
299 : !> \param i_offset block offset of atom_i's first RI sgf in int_3c dim 3 (default 0)
300 : !> \param screened .TRUE. if the kind-radius screens killed the whole triple (int_3c untouched)
301 : ! **************************************************************************************************
302 12994 : SUBROUTINE build_3c_integral_block_ctx(int_3c, ctx, ws, atom_j, atom_k, atom_i, &
303 : cell_j, cell_k, cell_i, &
304 : j_offset, k_offset, i_offset, screened)
305 :
306 : REAL(KIND=dp), DIMENSION(:, :, :) :: int_3c
307 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
308 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
309 : INTEGER, INTENT(IN) :: atom_j, atom_k, atom_i
310 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_j, cell_k, cell_i
311 : INTEGER, INTENT(IN), OPTIONAL :: j_offset, k_offset, i_offset
312 : LOGICAL, INTENT(OUT), OPTIONAL :: screened
313 :
314 : INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
315 : block_start_k, ikind, iset, jkind, jset, kkind, kset, my_i_offset, my_j_offset, &
316 : my_k_offset, ncoi, ncoj, ncok, nseti, nsetj, nsetk, sgfi, sgfj, sgfk
317 : INTEGER, DIMENSION(3) :: my_cell_i, my_cell_j, my_cell_k
318 12994 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
319 12994 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
320 12994 : nsgfj, nsgfk
321 12994 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
322 : REAL(KIND=dp) :: dij, dik, djk, kind_radius_i, &
323 : kind_radius_j, kind_radius_k, sijk_ext
324 12994 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sijk, sijk_contr
325 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
326 12994 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
327 12994 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, zeti, zetj, zetk
328 :
329 116 : IF (PRESENT(screened)) screened = .FALSE.
330 :
331 12994 : my_cell_i(1:3) = 0
332 12994 : IF (PRESENT(cell_i)) my_cell_i(1:3) = cell_i(1:3)
333 12994 : my_cell_j(1:3) = 0
334 12994 : IF (PRESENT(cell_j)) my_cell_j(1:3) = cell_j(1:3)
335 12994 : my_cell_k(1:3) = 0
336 12994 : IF (PRESENT(cell_k)) my_cell_k(1:3) = cell_k(1:3)
337 12994 : my_i_offset = 0
338 12994 : IF (PRESENT(i_offset)) my_i_offset = i_offset
339 12994 : my_j_offset = 0
340 12994 : IF (PRESENT(j_offset)) my_j_offset = j_offset
341 12994 : my_k_offset = 0
342 12994 : IF (PRESENT(k_offset)) my_k_offset = k_offset
343 :
344 259880 : ri = pbc(ctx%particle_set(atom_i)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_i, dp))
345 259880 : rj = pbc(ctx%particle_set(atom_j)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_j, dp))
346 259880 : rk = pbc(ctx%particle_set(atom_k)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_k, dp))
347 :
348 51976 : rjk(1:3) = rk(1:3) - rj(1:3)
349 51976 : rij(1:3) = rj(1:3) - ri(1:3)
350 51976 : rik(1:3) = rk(1:3) - ri(1:3)
351 :
352 51976 : djk = NORM2(rjk)
353 51976 : dij = NORM2(rij)
354 51976 : dik = NORM2(rik)
355 :
356 12994 : ikind = ctx%kind_of(atom_i)
357 12994 : jkind = ctx%kind_of(atom_j)
358 12994 : kkind = ctx%kind_of(atom_k)
359 :
360 : CALL get_gto_basis_set(ctx%basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, &
361 : lmax=lmax_i, lmin=lmin_i, npgf=npgfi, nset=nseti, &
362 : nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
363 12994 : zet=zeti, kind_radius=kind_radius_i)
364 : CALL get_gto_basis_set(ctx%basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, &
365 : lmax=lmax_j, lmin=lmin_j, npgf=npgfj, nset=nsetj, &
366 : nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
367 12994 : zet=zetj, kind_radius=kind_radius_j)
368 : CALL get_gto_basis_set(ctx%basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, &
369 : lmax=lmax_k, lmin=lmin_k, npgf=npgfk, nset=nsetk, &
370 : nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
371 12994 : zet=zetk, kind_radius=kind_radius_k)
372 :
373 : IF (kind_radius_j + kind_radius_i + ctx%dr_ij < dij .OR. &
374 12994 : kind_radius_j + kind_radius_k + ctx%dr_jk < djk .OR. &
375 : kind_radius_k + kind_radius_i + ctx%dr_ik < dik) THEN
376 0 : IF (PRESENT(screened)) screened = .TRUE.
377 0 : RETURN
378 : END IF
379 :
380 30595 : DO iset = 1, nseti
381 68170 : DO jset = 1, nsetj
382 37575 : IF (set_radius_j(jset) + set_radius_i(iset) + ctx%dr_ij < dij) CYCLE
383 124661 : DO kset = 1, nsetk
384 73229 : IF (set_radius_j(jset) + set_radius_k(kset) + ctx%dr_jk < djk) CYCLE
385 65129 : IF (set_radius_k(kset) + set_radius_i(iset) + ctx%dr_ik < dik) CYCLE
386 :
387 59587 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
388 59587 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
389 59587 : ncok = npgfk(kset)*ncoset(lmax_k(kset))
390 :
391 59587 : sgfi = first_sgf_i(1, iset)
392 59587 : sgfj = first_sgf_j(1, jset)
393 59587 : sgfk = first_sgf_k(1, kset)
394 :
395 59587 : IF (ncoj*ncok*ncoi <= 0) CYCLE
396 297935 : ALLOCATE (sijk(ncoj, ncok, ncoi))
397 59587 : sijk(:, :, :) = 0.0_dp
398 :
399 : CALL eri_3center(sijk, &
400 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
401 : rpgf_j(:, jset), rj, &
402 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), &
403 : rpgf_k(:, kset), rk, &
404 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
405 : rpgf_i(:, iset), ri, &
406 : djk, dij, dik, ws%lib, ctx%potential_parameter, &
407 59587 : int_abc_ext=sijk_ext)
408 :
409 297935 : ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
410 : CALL abc_contract_xsmm(sijk_contr, sijk, ctx%tspj(jset, jkind)%array, &
411 : ctx%spk(kset, kkind)%array, ctx%spi(iset, ikind)%array, &
412 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
413 59587 : nsgfi(iset), ws%cpp_buffer, ws%ccp_buffer)
414 59587 : DEALLOCATE (sijk)
415 :
416 59587 : block_start_j = sgfj + my_j_offset
417 59587 : block_end_j = sgfj + nsgfj(jset) - 1 + my_j_offset
418 59587 : block_start_k = sgfk + my_k_offset
419 59587 : block_end_k = sgfk + nsgfk(kset) - 1 + my_k_offset
420 59587 : block_start_i = sgfi + my_i_offset
421 59587 : block_end_i = sgfi + nsgfi(iset) - 1 + my_i_offset
422 :
423 : int_3c(block_start_j:block_end_j, &
424 : block_start_k:block_end_k, &
425 : block_start_i:block_end_i) = &
426 : int_3c(block_start_j:block_end_j, &
427 : block_start_k:block_end_k, &
428 : block_start_i:block_end_i) + &
429 722840 : sijk_contr(:, :, :)
430 110804 : DEALLOCATE (sijk_contr)
431 :
432 : END DO
433 : END DO
434 : END DO
435 :
436 38982 : END SUBROUTINE build_3c_integral_block_ctx
437 :
438 : ! **************************************************************************************************
439 : !> \brief ...
440 : !> \param int_3c ...
441 : !> \param qs_env ...
442 : !> \param potential_parameter ...
443 : !> \param basis_j ...
444 : !> \param basis_k ...
445 : !> \param basis_i ...
446 : !> \param cell_j ...
447 : !> \param cell_k ...
448 : !> \param cell_i ...
449 : !> \param atom_j ...
450 : !> \param atom_k ...
451 : !> \param atom_i ...
452 : !> \param j_bf_start_from_atom ...
453 : !> \param k_bf_start_from_atom ...
454 : !> \param i_bf_start_from_atom ...
455 : ! **************************************************************************************************
456 38550 : SUBROUTINE build_3c_integral_block(int_3c, qs_env, potential_parameter, &
457 12850 : basis_j, basis_k, basis_i, &
458 : cell_j, cell_k, cell_i, atom_j, atom_k, atom_i, &
459 12850 : j_bf_start_from_atom, k_bf_start_from_atom, &
460 12850 : i_bf_start_from_atom)
461 :
462 : REAL(KIND=dp), DIMENSION(:, :, :) :: int_3c
463 : TYPE(qs_environment_type), POINTER :: qs_env
464 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
465 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_j, basis_k, basis_i
466 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_j, cell_k, cell_i
467 : INTEGER, INTENT(IN), OPTIONAL :: atom_j, atom_k, atom_i
468 : INTEGER, DIMENSION(:), OPTIONAL :: j_bf_start_from_atom, &
469 : k_bf_start_from_atom, &
470 : i_bf_start_from_atom
471 :
472 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integral_block'
473 :
474 : INTEGER :: at_i, at_j, at_k, handle, my_i_offset, &
475 : my_j_offset, my_k_offset
476 167050 : TYPE(gw_3c_ctx_type) :: ctx
477 12850 : TYPE(gw_3c_ws_type) :: ws
478 :
479 12850 : CALL timeset(routineN, handle)
480 :
481 12850 : CALL gw_3c_ctx_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i)
482 12850 : CALL gw_3c_ws_create(ws, ctx)
483 :
484 680727 : int_3c(:, :, :) = 0.0_dp
485 :
486 : ! loop over all RI atoms
487 51087 : DO at_i = 1, ctx%natom
488 : ! loop over all AO atoms
489 165172 : DO at_j = 1, ctx%natom
490 : ! loop over all AO atoms
491 493325 : DO at_k = 1, ctx%natom
492 :
493 341003 : IF (PRESENT(atom_i)) THEN
494 340779 : IF (at_i /= atom_i) CYCLE
495 : END IF
496 114197 : IF (PRESENT(atom_j)) THEN
497 114197 : IF (at_j /= atom_j) CYCLE
498 : END IF
499 38293 : IF (PRESENT(atom_k)) THEN
500 38293 : IF (at_k /= atom_k) CYCLE
501 : END IF
502 :
503 12878 : IF (PRESENT(atom_j)) THEN
504 12878 : my_j_offset = 0
505 : ELSE
506 0 : CPASSERT(PRESENT(j_bf_start_from_atom))
507 0 : my_j_offset = j_bf_start_from_atom(at_j) - 1
508 : END IF
509 12878 : IF (PRESENT(atom_k)) THEN
510 12878 : my_k_offset = 0
511 : ELSE
512 0 : CPASSERT(PRESENT(k_bf_start_from_atom))
513 0 : my_k_offset = k_bf_start_from_atom(at_k) - 1
514 : END IF
515 12878 : IF (PRESENT(atom_i)) THEN
516 12822 : my_i_offset = 0
517 : ELSE
518 56 : CPASSERT(PRESENT(i_bf_start_from_atom))
519 56 : my_i_offset = i_bf_start_from_atom(at_i) - 1
520 : END IF
521 :
522 : CALL build_3c_integral_block_ctx(int_3c, ctx, ws, at_j, at_k, at_i, &
523 : cell_j=cell_j, cell_k=cell_k, cell_i=cell_i, &
524 : j_offset=my_j_offset, k_offset=my_k_offset, &
525 455088 : i_offset=my_i_offset)
526 :
527 : END DO ! atom_k (AO)
528 : END DO ! atom_j (AO)
529 : END DO ! atom_i (RI)
530 :
531 12850 : CALL gw_3c_ws_release(ws)
532 12850 : CALL gw_3c_ctx_release(ctx)
533 :
534 12850 : CALL timestop(handle)
535 :
536 12850 : END SUBROUTINE build_3c_integral_block
537 :
538 0 : END MODULE gw_integrals
539 :
|
