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_utils_compute_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 post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
42 : USE qs_environment_types, ONLY: get_qs_env,&
43 : qs_environment_type
44 : USE t_c_g0, ONLY: get_lmax_init,&
45 : init
46 :
47 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
48 : #include "./base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 :
52 : PRIVATE
53 :
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils_compute_integrals'
55 :
56 : PUBLIC :: build_3c_integral_block, build_3c_integral_block_ctx, &
57 : gw_3c_ctx_type, gw_3c_ctx_create, gw_3c_ctx_release, &
58 : gw_3c_ws_type, gw_3c_ws_create, gw_3c_ws_release
59 :
60 : ! **************************************************************************************************
61 : !> \brief Shared read-only context for repeated 3-center integral block builds: screening
62 : !> parameters, basis maxima, contracted sphi tables, and the one-time gamma /
63 : !> truncated-Coulomb table initializations. Create and release OUTSIDE any OMP parallel
64 : !> region; creation is MPI-collective when the potential is truncated.
65 : ! **************************************************************************************************
66 : TYPE gw_3c_ctx_type
67 : TYPE(libint_potential_type) :: potential_parameter = libint_potential_type()
68 : INTEGER :: op_ij = do_potential_id, &
69 : op_jk = do_potential_id
70 : REAL(KIND=dp) :: dr_ij = 0.0_dp, dr_jk = 0.0_dp, &
71 : dr_ik = 0.0_dp
72 : INTEGER :: maxli = 0, maxlj = 0, maxlk = 0, &
73 : max_am = 0, m_max = 0
74 : INTEGER :: max_ncoi = 0, max_ncoj = 0, max_ncok = 0
75 : INTEGER :: max_nsgfi = 0, max_nsgfj = 0, &
76 : max_nsgfk = 0, max_nset = 0, natom = 0
77 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi => NULL(), tspj => NULL(), &
78 : spk => NULL()
79 : TYPE(gto_basis_set_p_type), DIMENSION(:), ALLOCATABLE :: basis_i, basis_j, basis_k
80 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
81 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set => NULL()
82 : TYPE(cell_type), POINTER :: cell => NULL()
83 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat = 0.0_dp
84 : END TYPE gw_3c_ctx_type
85 :
86 : ! **************************************************************************************************
87 : !> \brief Per-thread workspace for 3-center integral block builds: libint object + contraction
88 : !> buffers. Each thread creates its own (inside the parallel region is fine).
89 : ! **************************************************************************************************
90 : TYPE gw_3c_ws_type
91 : TYPE(cp_libint_t), ALLOCATABLE :: lib
92 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpp_buffer, ccp_buffer
93 : END TYPE gw_3c_ws_type
94 :
95 : CONTAINS
96 :
97 : ! **************************************************************************************************
98 : !> \brief Build the shared 3c-integral context from the band-structure environment and explicitly
99 : !> supplied potential and basis sets.
100 : !> \param ctx ...
101 : !> \param bs_env ...
102 : !> \param potential_parameter ...
103 : !> \param basis_j ...
104 : !> \param basis_k ...
105 : !> \param basis_i ...
106 : ! **************************************************************************************************
107 780 : SUBROUTINE gw_3c_ctx_create(ctx, bs_env, potential_parameter, basis_j, basis_k, basis_i)
108 : TYPE(gw_3c_ctx_type), INTENT(OUT) :: ctx
109 : TYPE(post_scf_bandstructure_type), POINTER :: bs_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 60 : CPASSERT(ASSOCIATED(bs_env%ri_rs%atomic_kind_set))
114 60 : CPASSERT(ASSOCIATED(bs_env%ri_rs%cell))
115 60 : CPASSERT(ASSOCIATED(bs_env%para_env))
116 60 : CPASSERT(ASSOCIATED(bs_env%ri_rs%particle_set))
117 :
118 : CALL gw_3c_ctx_create_core(ctx, potential_parameter, basis_j, basis_k, basis_i, &
119 : bs_env%ri_rs%atomic_kind_set, bs_env%ri_rs%cell, &
120 60 : bs_env%n_atom, bs_env%para_env, bs_env%ri_rs%particle_set)
121 60 : END SUBROUTINE gw_3c_ctx_create
122 :
123 : ! **************************************************************************************************
124 : !> \brief Build screening and contraction data shared by repeated 3c-integral block evaluations.
125 : !> \param ctx ...
126 : !> \param potential_parameter ...
127 : !> \param basis_j ...
128 : !> \param basis_k ...
129 : !> \param basis_i ...
130 : !> \param atomic_kind_set ...
131 : !> \param cell ...
132 : !> \param natom ...
133 : !> \param para_env ...
134 : !> \param particle_set ...
135 : ! **************************************************************************************************
136 181748 : SUBROUTINE gw_3c_ctx_create_core(ctx, potential_parameter, basis_j, basis_k, basis_i, &
137 : atomic_kind_set, cell, natom, para_env, particle_set)
138 : TYPE(gw_3c_ctx_type), INTENT(OUT) :: ctx
139 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
140 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_j, basis_k, basis_i
141 : TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), &
142 : POINTER :: atomic_kind_set
143 : TYPE(cell_type), INTENT(IN), POINTER :: cell
144 : INTEGER, INTENT(IN) :: natom
145 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
146 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
147 : POINTER :: particle_set
148 :
149 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_3c_ctx_create_core'
150 :
151 : INTEGER :: egfi, handle, ibasis, ilist, imax, iset, &
152 : jset, kset, l, nbasis, ncoi, npgf_l, &
153 : sgfi, unit_id
154 12982 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmin_i, npgfi, npgfj, npgfk, &
155 12982 : nsgfi, nsgfj, nsgfk
156 : TYPE(gto_basis_set_type), POINTER :: basis_set
157 :
158 12982 : CALL timeset(routineN, handle)
159 :
160 12982 : ctx%potential_parameter = potential_parameter
161 12982 : ctx%op_ij = potential_parameter%potential_type
162 12982 : ctx%op_jk = do_potential_id
163 :
164 12982 : IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_ij == do_potential_short) THEN
165 12982 : ctx%dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
166 12982 : ctx%dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
167 0 : ELSE IF (ctx%op_ij == do_potential_coulomb) THEN
168 0 : ctx%dr_ij = 1000000.0_dp
169 0 : ctx%dr_ik = 1000000.0_dp
170 : END IF
171 :
172 12982 : ctx%cell => cell
173 12982 : ctx%natom = natom
174 12982 : ctx%particle_set => particle_set
175 12982 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=ctx%kind_of)
176 12982 : CALL get_cell(cell=ctx%cell, h=ctx%hmat)
177 :
178 51802 : ctx%basis_i = basis_i
179 51802 : ctx%basis_j = basis_j
180 51802 : ctx%basis_k = basis_k
181 :
182 : ! max l per basis for libint; max nset/nco/nsgf for the LIBXSMM contraction buffers
183 12982 : nbasis = SIZE(basis_i)
184 38820 : DO ibasis = 1, nbasis
185 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
186 25838 : lmax=lmax_i, lmin=lmin_i, nset=iset, nsgf_set=nsgfi, npgf=npgfi)
187 25838 : ctx%maxli = MAX(ctx%maxli, imax)
188 25838 : ctx%max_nset = MAX(ctx%max_nset, iset)
189 66002 : ctx%max_nsgfi = MAX(ctx%max_nsgfi, MAXVAL(nsgfi))
190 66002 : ctx%max_ncoi = MAX(ctx%max_ncoi, MAXVAL(npgfi)*ncoset(ctx%maxli))
191 103716 : DO l = 0, imax
192 39058 : npgf_l = 0
193 109424 : DO jset = 1, iset
194 109424 : IF (lmin_i(jset) == l .AND. lmax_i(jset) == l) npgf_l = npgf_l + npgfi(jset)
195 : END DO
196 64896 : ctx%max_ncoi = MAX(ctx%max_ncoi, npgf_l*ncoset(l))
197 : END DO
198 : END DO
199 38820 : DO ibasis = 1, nbasis
200 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
201 25838 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
202 25838 : ctx%maxlj = MAX(ctx%maxlj, imax)
203 25838 : ctx%max_nset = MAX(ctx%max_nset, jset)
204 77640 : ctx%max_nsgfj = MAX(ctx%max_nsgfj, MAXVAL(nsgfj))
205 116460 : ctx%max_ncoj = MAX(ctx%max_ncoj, MAXVAL(npgfj)*ncoset(ctx%maxlj))
206 : END DO
207 38820 : DO ibasis = 1, nbasis
208 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
209 25838 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
210 25838 : ctx%maxlk = MAX(ctx%maxlk, imax)
211 25838 : ctx%max_nset = MAX(ctx%max_nset, kset)
212 77640 : ctx%max_nsgfk = MAX(ctx%max_nsgfk, MAXVAL(nsgfk))
213 116460 : ctx%max_ncok = MAX(ctx%max_ncok, MAXVAL(npgfk)*ncoset(ctx%maxlk))
214 : END DO
215 12982 : ctx%m_max = ctx%maxli + ctx%maxlj + ctx%maxlk
216 12982 : ctx%max_am = MAX(ctx%maxli, ctx%maxlj, ctx%maxlk)
217 :
218 : ! contiguous (and for j transposed) sphi copies, shared read-only across threads
219 : ALLOCATE (ctx%spi(ctx%max_nset, nbasis), ctx%tspj(ctx%max_nset, nbasis), &
220 443824 : ctx%spk(ctx%max_nset, nbasis))
221 38820 : DO ibasis = 1, nbasis
222 117650 : DO iset = 1, ctx%max_nset
223 78830 : NULLIFY (ctx%spi(iset, ibasis)%array)
224 78830 : NULLIFY (ctx%tspj(iset, ibasis)%array)
225 104668 : NULLIFY (ctx%spk(iset, ibasis)%array)
226 : END DO
227 : END DO
228 51928 : DO ilist = 1, 3
229 129442 : DO ibasis = 1, nbasis
230 77514 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
231 77514 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
232 77514 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
233 260228 : DO iset = 1, basis_set%nset
234 143768 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
235 143768 : sgfi = basis_set%first_sgf(1, iset)
236 143768 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
237 221282 : IF (ilist == 1) THEN
238 160656 : ALLOCATE (ctx%spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
239 346322 : ctx%spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
240 103604 : ELSE IF (ilist == 2) THEN
241 207208 : ALLOCATE (ctx%tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
242 1646902 : ctx%tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
243 : ELSE
244 207208 : ALLOCATE (ctx%spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
245 1402574 : ctx%spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
246 : END IF
247 : END DO
248 : END DO
249 : END DO
250 :
251 : ! one-time table inits; the truncated-Coulomb init reads a file + bcasts => MPI-collective,
252 : ! must happen here and never inside the per-block path or an OMP region
253 12982 : IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_jk == do_potential_truncated) THEN
254 12982 : IF (ctx%m_max > get_lmax_init()) THEN
255 0 : IF (para_env%mepos == 0) THEN
256 0 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
257 : END IF
258 0 : CALL init(ctx%m_max, unit_id, para_env%mepos, para_env)
259 0 : IF (para_env%mepos == 0) THEN
260 0 : CALL close_file(unit_id)
261 : END IF
262 : END IF
263 : END IF
264 12982 : CALL init_md_ftable(nmax=ctx%m_max)
265 :
266 12982 : CALL timestop(handle)
267 :
268 25964 : END SUBROUTINE gw_3c_ctx_create_core
269 :
270 : ! **************************************************************************************************
271 : !> \brief Releases the shared 3c-integral context.
272 : !> \param ctx ...
273 : ! **************************************************************************************************
274 12982 : SUBROUTINE gw_3c_ctx_release(ctx)
275 :
276 : TYPE(gw_3c_ctx_type), INTENT(INOUT) :: ctx
277 :
278 : INTEGER :: ibasis, iset
279 :
280 52900 : DO iset = 1, SIZE(ctx%spi, 1)
281 131730 : DO ibasis = 1, SIZE(ctx%spi, 2)
282 78830 : IF (ASSOCIATED(ctx%spi(iset, ibasis)%array)) DEALLOCATE (ctx%spi(iset, ibasis)%array)
283 78830 : IF (ASSOCIATED(ctx%tspj(iset, ibasis)%array)) DEALLOCATE (ctx%tspj(iset, ibasis)%array)
284 118748 : IF (ASSOCIATED(ctx%spk(iset, ibasis)%array)) DEALLOCATE (ctx%spk(iset, ibasis)%array)
285 : END DO
286 : END DO
287 12982 : DEALLOCATE (ctx%spi, ctx%tspj, ctx%spk)
288 12982 : NULLIFY (ctx%spi, ctx%tspj, ctx%spk, ctx%particle_set, ctx%cell)
289 12982 : IF (ALLOCATED(ctx%kind_of)) DEALLOCATE (ctx%kind_of)
290 12982 : IF (ALLOCATED(ctx%basis_i)) DEALLOCATE (ctx%basis_i)
291 12982 : IF (ALLOCATED(ctx%basis_j)) DEALLOCATE (ctx%basis_j)
292 12982 : IF (ALLOCATED(ctx%basis_k)) DEALLOCATE (ctx%basis_k)
293 :
294 12982 : END SUBROUTINE gw_3c_ctx_release
295 :
296 : ! **************************************************************************************************
297 : !> \brief Creates a per-thread 3c workspace: libint object + LIBXSMM contraction buffers.
298 : !> \param ws ...
299 : !> \param ctx ...
300 : ! **************************************************************************************************
301 12992 : SUBROUTINE gw_3c_ws_create(ws, ctx)
302 :
303 : TYPE(gw_3c_ws_type), INTENT(OUT) :: ws
304 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
305 :
306 12992 : ALLOCATE (ws%lib)
307 12992 : CALL cp_libint_init_3eri(ws%lib, ctx%max_am)
308 12992 : CALL cp_libint_set_contrdepth(ws%lib, 1)
309 0 : ALLOCATE (ws%cpp_buffer(ctx%max_nsgfj*ctx%max_ncok), &
310 64960 : ws%ccp_buffer(ctx%max_nsgfj*ctx%max_nsgfk*ctx%max_ncoi))
311 :
312 12992 : END SUBROUTINE gw_3c_ws_create
313 :
314 : ! **************************************************************************************************
315 : !> \brief Releases a per-thread 3c workspace.
316 : !> \param ws ...
317 : ! **************************************************************************************************
318 12992 : SUBROUTINE gw_3c_ws_release(ws)
319 :
320 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
321 :
322 12992 : CALL cp_libint_cleanup_3eri(ws%lib)
323 12992 : DEALLOCATE (ws%lib)
324 12992 : DEALLOCATE (ws%cpp_buffer, ws%ccp_buffer)
325 :
326 12992 : END SUBROUTINE gw_3c_ws_release
327 :
328 : ! **************************************************************************************************
329 : !> \brief Computes the 3c integral block (mu(atom_j) nu(atom_k) | P(atom_i)) for ONE atom triple,
330 : !> accumulating into the caller-zeroed int_3c at the given block offsets.
331 : !> Thread-safe: reads the frozen ctx + read-only module tables, mutates only its arguments
332 : !> and the per-thread ws.
333 : !> \param int_3c pre-zeroed target; the triple's contribution is accumulated in place
334 : !> \param ctx shared context from gw_3c_ctx_create
335 : !> \param ws per-thread workspace from gw_3c_ws_create
336 : !> \param atom_j ...
337 : !> \param atom_k ...
338 : !> \param atom_i ...
339 : !> \param cell_j ...
340 : !> \param cell_k ...
341 : !> \param cell_i ...
342 : !> \param j_offset block offset of atom_j's first sgf in int_3c dim 1 (default 0)
343 : !> \param k_offset block offset of atom_k's first sgf in int_3c dim 2 (default 0)
344 : !> \param i_offset block offset of atom_i's first RI sgf in int_3c dim 3 (default 0)
345 : !> \param screened .TRUE. if the kind-radius screens killed the whole triple (int_3c untouched)
346 : ! **************************************************************************************************
347 13673 : SUBROUTINE build_3c_integral_block_ctx(int_3c, ctx, ws, atom_j, atom_k, atom_i, &
348 : cell_j, cell_k, cell_i, &
349 : j_offset, k_offset, i_offset, screened)
350 :
351 : REAL(KIND=dp), DIMENSION(:, :, :) :: int_3c
352 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
353 : TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
354 : INTEGER, INTENT(IN) :: atom_j, atom_k, atom_i
355 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_j, cell_k, cell_i
356 : INTEGER, INTENT(IN), OPTIONAL :: j_offset, k_offset, i_offset
357 : LOGICAL, INTENT(OUT), OPTIONAL :: screened
358 :
359 : INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
360 : block_start_k, ikind, iset, jkind, jset, kkind, kset, my_i_offset, my_j_offset, &
361 : my_k_offset, ncoi, ncoj, ncok, nseti, nsetj, nsetk, sgfi, sgfj, sgfk
362 : INTEGER, DIMENSION(3) :: my_cell_i, my_cell_j, my_cell_k
363 13673 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
364 13673 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
365 13673 : nsgfj, nsgfk
366 13673 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
367 : REAL(KIND=dp) :: dij, dik, djk, kind_radius_i, &
368 : kind_radius_j, kind_radius_k, sijk_ext
369 13673 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sijk, sijk_contr
370 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
371 13673 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
372 13673 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, zeti, zetj, zetk
373 :
374 651 : IF (PRESENT(screened)) screened = .FALSE.
375 :
376 13673 : my_cell_i(1:3) = 0
377 13673 : IF (PRESENT(cell_i)) my_cell_i(1:3) = cell_i(1:3)
378 13673 : my_cell_j(1:3) = 0
379 13673 : IF (PRESENT(cell_j)) my_cell_j(1:3) = cell_j(1:3)
380 13673 : my_cell_k(1:3) = 0
381 13673 : IF (PRESENT(cell_k)) my_cell_k(1:3) = cell_k(1:3)
382 13673 : my_i_offset = 0
383 13673 : IF (PRESENT(i_offset)) my_i_offset = i_offset
384 13673 : my_j_offset = 0
385 13673 : IF (PRESENT(j_offset)) my_j_offset = j_offset
386 13673 : my_k_offset = 0
387 13673 : IF (PRESENT(k_offset)) my_k_offset = k_offset
388 :
389 273460 : ri = pbc(ctx%particle_set(atom_i)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_i, dp))
390 273460 : rj = pbc(ctx%particle_set(atom_j)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_j, dp))
391 273460 : rk = pbc(ctx%particle_set(atom_k)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_k, dp))
392 :
393 54692 : rjk(1:3) = rk(1:3) - rj(1:3)
394 54692 : rij(1:3) = rj(1:3) - ri(1:3)
395 54692 : rik(1:3) = rk(1:3) - ri(1:3)
396 :
397 54692 : djk = NORM2(rjk)
398 54692 : dij = NORM2(rij)
399 54692 : dik = NORM2(rik)
400 :
401 13673 : ikind = ctx%kind_of(atom_i)
402 13673 : jkind = ctx%kind_of(atom_j)
403 13673 : kkind = ctx%kind_of(atom_k)
404 :
405 : CALL get_gto_basis_set(ctx%basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, &
406 : lmax=lmax_i, lmin=lmin_i, npgf=npgfi, nset=nseti, &
407 : nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
408 13673 : zet=zeti, kind_radius=kind_radius_i)
409 : CALL get_gto_basis_set(ctx%basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, &
410 : lmax=lmax_j, lmin=lmin_j, npgf=npgfj, nset=nsetj, &
411 : nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
412 13673 : zet=zetj, kind_radius=kind_radius_j)
413 : CALL get_gto_basis_set(ctx%basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, &
414 : lmax=lmax_k, lmin=lmin_k, npgf=npgfk, nset=nsetk, &
415 : nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
416 13673 : zet=zetk, kind_radius=kind_radius_k)
417 :
418 : IF (kind_radius_j + kind_radius_i + ctx%dr_ij < dij .OR. &
419 13673 : kind_radius_j + kind_radius_k + ctx%dr_jk < djk .OR. &
420 : kind_radius_k + kind_radius_i + ctx%dr_ik < dik) THEN
421 0 : IF (PRESENT(screened)) screened = .TRUE.
422 0 : RETURN
423 : END IF
424 :
425 37803 : DO iset = 1, nseti
426 89242 : DO jset = 1, nsetj
427 51439 : IF (set_radius_j(jset) + set_radius_i(iset) + ctx%dr_ij < dij) CYCLE
428 183684 : DO kset = 1, nsetk
429 111859 : IF (set_radius_j(jset) + set_radius_k(kset) + ctx%dr_jk < djk) CYCLE
430 103759 : IF (set_radius_k(kset) + set_radius_i(iset) + ctx%dr_ik < dik) CYCLE
431 :
432 98217 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
433 98217 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
434 98217 : ncok = npgfk(kset)*ncoset(lmax_k(kset))
435 :
436 98217 : sgfi = first_sgf_i(1, iset)
437 98217 : sgfj = first_sgf_j(1, jset)
438 98217 : sgfk = first_sgf_k(1, kset)
439 :
440 98217 : IF (ncoj*ncok*ncoi <= 0) CYCLE
441 491085 : ALLOCATE (sijk(ncoj, ncok, ncoi))
442 98217 : sijk(:, :, :) = 0.0_dp
443 :
444 : CALL eri_3center(sijk, &
445 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
446 : rpgf_j(:, jset), rj, &
447 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), &
448 : rpgf_k(:, kset), rk, &
449 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
450 : rpgf_i(:, iset), ri, &
451 : djk, dij, dik, ws%lib, ctx%potential_parameter, &
452 98217 : int_abc_ext=sijk_ext)
453 :
454 491085 : ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
455 : CALL abc_contract_xsmm(sijk_contr, sijk, ctx%tspj(jset, jkind)%array, &
456 : ctx%spk(kset, kkind)%array, ctx%spi(iset, ikind)%array, &
457 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
458 98217 : nsgfi(iset), ws%cpp_buffer, ws%ccp_buffer)
459 98217 : DEALLOCATE (sijk)
460 :
461 98217 : block_start_j = sgfj + my_j_offset
462 98217 : block_end_j = sgfj + nsgfj(jset) - 1 + my_j_offset
463 98217 : block_start_k = sgfk + my_k_offset
464 98217 : block_end_k = sgfk + nsgfk(kset) - 1 + my_k_offset
465 98217 : block_start_i = sgfi + my_i_offset
466 98217 : block_end_i = sgfi + nsgfi(iset) - 1 + my_i_offset
467 :
468 : int_3c(block_start_j:block_end_j, &
469 : block_start_k:block_end_k, &
470 : block_start_i:block_end_i) = &
471 : int_3c(block_start_j:block_end_j, &
472 : block_start_k:block_end_k, &
473 : block_start_i:block_end_i) + &
474 1308301 : sijk_contr(:, :, :)
475 163298 : DEALLOCATE (sijk_contr)
476 :
477 : END DO
478 : END DO
479 : END DO
480 :
481 41019 : END SUBROUTINE build_3c_integral_block_ctx
482 :
483 : ! **************************************************************************************************
484 : !> \brief ...
485 : !> \param int_3c ...
486 : !> \param qs_env ...
487 : !> \param potential_parameter ...
488 : !> \param basis_j ...
489 : !> \param basis_k ...
490 : !> \param basis_i ...
491 : !> \param cell_j ...
492 : !> \param cell_k ...
493 : !> \param cell_i ...
494 : !> \param atom_j ...
495 : !> \param atom_k ...
496 : !> \param atom_i ...
497 : !> \param j_bf_start_from_atom ...
498 : !> \param k_bf_start_from_atom ...
499 : !> \param i_bf_start_from_atom ...
500 : ! **************************************************************************************************
501 38766 : SUBROUTINE build_3c_integral_block(int_3c, qs_env, potential_parameter, &
502 12922 : basis_j, basis_k, basis_i, &
503 : cell_j, cell_k, cell_i, atom_j, atom_k, atom_i, &
504 12922 : j_bf_start_from_atom, k_bf_start_from_atom, &
505 12922 : i_bf_start_from_atom)
506 :
507 : REAL(KIND=dp), DIMENSION(:, :, :) :: int_3c
508 : TYPE(qs_environment_type), POINTER :: qs_env
509 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
510 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_j, basis_k, basis_i
511 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_j, cell_k, cell_i
512 : INTEGER, INTENT(IN), OPTIONAL :: atom_j, atom_k, atom_i
513 : INTEGER, DIMENSION(:), OPTIONAL :: j_bf_start_from_atom, &
514 : k_bf_start_from_atom, &
515 : i_bf_start_from_atom
516 :
517 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integral_block'
518 :
519 : INTEGER :: at_i, at_j, at_k, handle, my_i_offset, &
520 : my_j_offset, my_k_offset, natom
521 12922 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
522 : TYPE(cell_type), POINTER :: cell
523 167986 : TYPE(gw_3c_ctx_type) :: ctx
524 12922 : TYPE(gw_3c_ws_type) :: ws
525 : TYPE(mp_para_env_type), POINTER :: para_env
526 12922 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
527 :
528 12922 : CALL timeset(routineN, handle)
529 :
530 12922 : NULLIFY (atomic_kind_set, cell, para_env, particle_set)
531 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=cell, natom=natom, &
532 12922 : para_env=para_env, particle_set=particle_set)
533 : CALL gw_3c_ctx_create_core(ctx, potential_parameter, basis_j, basis_k, basis_i, &
534 12922 : atomic_kind_set, cell, natom, para_env, particle_set)
535 12922 : CALL gw_3c_ws_create(ws, ctx)
536 :
537 823215 : int_3c(:, :, :) = 0.0_dp
538 :
539 : ! loop over all RI atoms
540 51303 : DO at_i = 1, ctx%natom
541 : ! loop over all AO atoms
542 165676 : DO at_j = 1, ctx%natom
543 : ! loop over all AO atoms
544 494333 : DO at_k = 1, ctx%natom
545 :
546 341579 : IF (PRESENT(atom_i)) THEN
547 340779 : IF (at_i /= atom_i) CYCLE
548 : END IF
549 114773 : IF (PRESENT(atom_j)) THEN
550 114773 : IF (at_j /= atom_j) CYCLE
551 : END IF
552 38581 : IF (PRESENT(atom_k)) THEN
553 38581 : IF (at_k /= atom_k) CYCLE
554 : END IF
555 :
556 13022 : IF (PRESENT(atom_j)) THEN
557 13022 : my_j_offset = 0
558 : ELSE
559 0 : CPASSERT(PRESENT(j_bf_start_from_atom))
560 0 : my_j_offset = j_bf_start_from_atom(at_j) - 1
561 : END IF
562 13022 : IF (PRESENT(atom_k)) THEN
563 13022 : my_k_offset = 0
564 : ELSE
565 0 : CPASSERT(PRESENT(k_bf_start_from_atom))
566 0 : my_k_offset = k_bf_start_from_atom(at_k) - 1
567 : END IF
568 13022 : IF (PRESENT(atom_i)) THEN
569 12822 : my_i_offset = 0
570 : ELSE
571 200 : CPASSERT(PRESENT(i_bf_start_from_atom))
572 200 : my_i_offset = i_bf_start_from_atom(at_i) - 1
573 : END IF
574 :
575 : CALL build_3c_integral_block_ctx(int_3c, ctx, ws, at_j, at_k, at_i, &
576 : cell_j=cell_j, cell_k=cell_k, cell_i=cell_i, &
577 : j_offset=my_j_offset, k_offset=my_k_offset, &
578 455952 : i_offset=my_i_offset)
579 :
580 : END DO ! atom_k (AO)
581 : END DO ! atom_j (AO)
582 : END DO ! atom_i (RI)
583 :
584 12922 : CALL gw_3c_ws_release(ws)
585 12922 : CALL gw_3c_ctx_release(ctx)
586 :
587 12922 : CALL timestop(handle)
588 :
589 38766 : END SUBROUTINE build_3c_integral_block
590 :
591 0 : END MODULE gw_utils_compute_integrals
|
