(git:1e0ab71)
Loading...
Searching...
No Matches
gw_integrals.F
Go to the documentation of this file.
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! **************************************************************************************************
11MODULE 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
94CONTAINS
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 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 INTEGER, DIMENSION(:), POINTER :: npgfi, npgfj, npgfk, nsgfi, nsgfj, nsgfk
118 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 CALL timeset(routinen, handle)
123
124 ctx%potential_parameter = potential_parameter
125 ctx%op_ij = potential_parameter%potential_type
126 ctx%op_jk = do_potential_id
127
128 IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_ij == do_potential_short) THEN
129 ctx%dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
130 ctx%dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
131 ELSEIF (ctx%op_ij == do_potential_coulomb) THEN
132 ctx%dr_ij = 1000000.0_dp
133 ctx%dr_ik = 1000000.0_dp
134 END IF
135
136 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 para_env=para_env, particle_set=ctx%particle_set)
139 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=ctx%kind_of)
140 CALL get_cell(cell=ctx%cell, h=ctx%hmat)
141
142 ctx%basis_i = basis_i
143 ctx%basis_j = basis_j
144 ctx%basis_k = basis_k
145
146 ! max l per basis for libint; max nset/nco/nsgf for the LIBXSMM contraction buffers
147 nbasis = SIZE(basis_i)
148 DO ibasis = 1, nbasis
149 CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
150 nset=iset, nsgf_set=nsgfi, npgf=npgfi)
151 ctx%maxli = max(ctx%maxli, imax)
152 ctx%max_nset = max(ctx%max_nset, iset)
153 ctx%max_nsgfi = max(ctx%max_nsgfi, maxval(nsgfi))
154 ctx%max_ncoi = max(ctx%max_ncoi, maxval(npgfi)*ncoset(ctx%maxli))
155 END DO
156 DO ibasis = 1, nbasis
157 CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
158 nset=jset, nsgf_set=nsgfj, npgf=npgfj)
159 ctx%maxlj = max(ctx%maxlj, imax)
160 ctx%max_nset = max(ctx%max_nset, jset)
161 ctx%max_nsgfj = max(ctx%max_nsgfj, maxval(nsgfj))
162 ctx%max_ncoj = max(ctx%max_ncoj, maxval(npgfj)*ncoset(ctx%maxlj))
163 END DO
164 DO ibasis = 1, nbasis
165 CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
166 nset=kset, nsgf_set=nsgfk, npgf=npgfk)
167 ctx%maxlk = max(ctx%maxlk, imax)
168 ctx%max_nset = max(ctx%max_nset, kset)
169 ctx%max_nsgfk = max(ctx%max_nsgfk, maxval(nsgfk))
170 ctx%max_ncok = max(ctx%max_ncok, maxval(npgfk)*ncoset(ctx%maxlk))
171 END DO
172 ctx%m_max = ctx%maxli + ctx%maxlj + ctx%maxlk
173 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 ctx%spk(ctx%max_nset, nbasis))
178 DO ibasis = 1, nbasis
179 DO iset = 1, ctx%max_nset
180 NULLIFY (ctx%spi(iset, ibasis)%array)
181 NULLIFY (ctx%tspj(iset, ibasis)%array)
182 NULLIFY (ctx%spk(iset, ibasis)%array)
183 END DO
184 END DO
185 DO ilist = 1, 3
186 DO ibasis = 1, nbasis
187 IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
188 IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
189 IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
190 DO iset = 1, basis_set%nset
191 ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
192 sgfi = basis_set%first_sgf(1, iset)
193 egfi = sgfi + basis_set%nsgf_set(iset) - 1
194 IF (ilist == 1) THEN
195 ALLOCATE (ctx%spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
196 ctx%spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
197 ELSE IF (ilist == 2) THEN
198 ALLOCATE (ctx%tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
199 ctx%tspj(iset, ibasis)%array(:, :) = transpose(basis_set%sphi(1:ncoi, sgfi:egfi))
200 ELSE
201 ALLOCATE (ctx%spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
202 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 IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_jk == do_potential_truncated) THEN
211 IF (ctx%m_max > get_lmax_init()) THEN
212 IF (para_env%mepos == 0) THEN
213 CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
214 END IF
215 CALL init(ctx%m_max, unit_id, para_env%mepos, para_env)
216 IF (para_env%mepos == 0) THEN
217 CALL close_file(unit_id)
218 END IF
219 END IF
220 END IF
221 CALL init_md_ftable(nmax=ctx%m_max)
222
223 CALL timestop(handle)
224
225 END SUBROUTINE gw_3c_ctx_create
226
227! **************************************************************************************************
228!> \brief Releases the shared 3c-integral context.
229!> \param ctx ...
230! **************************************************************************************************
231 SUBROUTINE gw_3c_ctx_release(ctx)
232
233 TYPE(gw_3c_ctx_type), INTENT(INOUT) :: ctx
234
235 INTEGER :: ibasis, iset
236
237 DO iset = 1, SIZE(ctx%spi, 1)
238 DO ibasis = 1, SIZE(ctx%spi, 2)
239 IF (ASSOCIATED(ctx%spi(iset, ibasis)%array)) DEALLOCATE (ctx%spi(iset, ibasis)%array)
240 IF (ASSOCIATED(ctx%tspj(iset, ibasis)%array)) DEALLOCATE (ctx%tspj(iset, ibasis)%array)
241 IF (ASSOCIATED(ctx%spk(iset, ibasis)%array)) DEALLOCATE (ctx%spk(iset, ibasis)%array)
242 END DO
243 END DO
244 DEALLOCATE (ctx%spi, ctx%tspj, ctx%spk)
245 NULLIFY (ctx%spi, ctx%tspj, ctx%spk, ctx%particle_set, ctx%cell)
246 IF (ALLOCATED(ctx%kind_of)) DEALLOCATE (ctx%kind_of)
247 IF (ALLOCATED(ctx%basis_i)) DEALLOCATE (ctx%basis_i)
248 IF (ALLOCATED(ctx%basis_j)) DEALLOCATE (ctx%basis_j)
249 IF (ALLOCATED(ctx%basis_k)) DEALLOCATE (ctx%basis_k)
250
251 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 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 CALL cp_libint_init_3eri(ws%lib, ctx%max_am)
264 CALL cp_libint_set_contrdepth(ws%lib, 1)
265 ALLOCATE (ws%cpp_buffer(ctx%max_nsgfj*ctx%max_ncok), &
266 ws%ccp_buffer(ctx%max_nsgfj*ctx%max_nsgfk*ctx%max_ncoi))
267
268 END SUBROUTINE gw_3c_ws_create
269
270! **************************************************************************************************
271!> \brief Releases a per-thread 3c workspace.
272!> \param ws ...
273! **************************************************************************************************
274 SUBROUTINE gw_3c_ws_release(ws)
275
276 TYPE(gw_3c_ws_type), INTENT(INOUT) :: ws
277
278 CALL cp_libint_cleanup_3eri(ws%lib)
279 DEALLOCATE (ws%cpp_buffer, ws%ccp_buffer)
280
281 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 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 INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
319 lmin_k, npgfi, npgfj, npgfk, nsgfi, &
320 nsgfj, nsgfk
321 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sijk, sijk_contr
325 REAL(kind=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
326 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
327 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, zeti, zetj, zetk
328
329 IF (PRESENT(screened)) screened = .false.
330
331 my_cell_i(1:3) = 0
332 IF (PRESENT(cell_i)) my_cell_i(1:3) = cell_i(1:3)
333 my_cell_j(1:3) = 0
334 IF (PRESENT(cell_j)) my_cell_j(1:3) = cell_j(1:3)
335 my_cell_k(1:3) = 0
336 IF (PRESENT(cell_k)) my_cell_k(1:3) = cell_k(1:3)
337 my_i_offset = 0
338 IF (PRESENT(i_offset)) my_i_offset = i_offset
339 my_j_offset = 0
340 IF (PRESENT(j_offset)) my_j_offset = j_offset
341 my_k_offset = 0
342 IF (PRESENT(k_offset)) my_k_offset = k_offset
343
344 ri = pbc(ctx%particle_set(atom_i)%r(1:3), ctx%cell) + matmul(ctx%hmat, real(my_cell_i, dp))
345 rj = pbc(ctx%particle_set(atom_j)%r(1:3), ctx%cell) + matmul(ctx%hmat, real(my_cell_j, dp))
346 rk = pbc(ctx%particle_set(atom_k)%r(1:3), ctx%cell) + matmul(ctx%hmat, real(my_cell_k, dp))
347
348 rjk(1:3) = rk(1:3) - rj(1:3)
349 rij(1:3) = rj(1:3) - ri(1:3)
350 rik(1:3) = rk(1:3) - ri(1:3)
351
352 djk = norm2(rjk)
353 dij = norm2(rij)
354 dik = norm2(rik)
355
356 ikind = ctx%kind_of(atom_i)
357 jkind = ctx%kind_of(atom_j)
358 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 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 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 zet=zetk, kind_radius=kind_radius_k)
372
373 IF (kind_radius_j + kind_radius_i + ctx%dr_ij < dij .OR. &
374 kind_radius_j + kind_radius_k + ctx%dr_jk < djk .OR. &
375 kind_radius_k + kind_radius_i + ctx%dr_ik < dik) THEN
376 IF (PRESENT(screened)) screened = .true.
377 RETURN
378 END IF
379
380 DO iset = 1, nseti
381 DO jset = 1, nsetj
382 IF (set_radius_j(jset) + set_radius_i(iset) + ctx%dr_ij < dij) cycle
383 DO kset = 1, nsetk
384 IF (set_radius_j(jset) + set_radius_k(kset) + ctx%dr_jk < djk) cycle
385 IF (set_radius_k(kset) + set_radius_i(iset) + ctx%dr_ik < dik) cycle
386
387 ncoi = npgfi(iset)*ncoset(lmax_i(iset))
388 ncoj = npgfj(jset)*ncoset(lmax_j(jset))
389 ncok = npgfk(kset)*ncoset(lmax_k(kset))
390
391 sgfi = first_sgf_i(1, iset)
392 sgfj = first_sgf_j(1, jset)
393 sgfk = first_sgf_k(1, kset)
394
395 IF (ncoj*ncok*ncoi <= 0) cycle
396 ALLOCATE (sijk(ncoj, ncok, ncoi))
397 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 int_abc_ext=sijk_ext)
408
409 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 nsgfi(iset), ws%cpp_buffer, ws%ccp_buffer)
414 DEALLOCATE (sijk)
415
416 block_start_j = sgfj + my_j_offset
417 block_end_j = sgfj + nsgfj(jset) - 1 + my_j_offset
418 block_start_k = sgfk + my_k_offset
419 block_end_k = sgfk + nsgfk(kset) - 1 + my_k_offset
420 block_start_i = sgfi + my_i_offset
421 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 sijk_contr(:, :, :)
430 DEALLOCATE (sijk_contr)
431
432 END DO
433 END DO
434 END DO
435
436 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 SUBROUTINE build_3c_integral_block(int_3c, qs_env, potential_parameter, &
457 basis_j, basis_k, basis_i, &
458 cell_j, cell_k, cell_i, atom_j, atom_k, atom_i, &
459 j_bf_start_from_atom, k_bf_start_from_atom, &
460 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 TYPE(gw_3c_ctx_type) :: ctx
477 TYPE(gw_3c_ws_type) :: ws
478
479 CALL timeset(routinen, handle)
480
481 CALL gw_3c_ctx_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i)
482 CALL gw_3c_ws_create(ws, ctx)
483
484 int_3c(:, :, :) = 0.0_dp
485
486 ! loop over all RI atoms
487 DO at_i = 1, ctx%natom
488 ! loop over all AO atoms
489 DO at_j = 1, ctx%natom
490 ! loop over all AO atoms
491 DO at_k = 1, ctx%natom
492
493 IF (PRESENT(atom_i)) THEN
494 IF (at_i /= atom_i) cycle
495 END IF
496 IF (PRESENT(atom_j)) THEN
497 IF (at_j /= atom_j) cycle
498 END IF
499 IF (PRESENT(atom_k)) THEN
500 IF (at_k /= atom_k) cycle
501 END IF
502
503 IF (PRESENT(atom_j)) THEN
504 my_j_offset = 0
505 ELSE
506 cpassert(PRESENT(j_bf_start_from_atom))
507 my_j_offset = j_bf_start_from_atom(at_j) - 1
508 END IF
509 IF (PRESENT(atom_k)) THEN
510 my_k_offset = 0
511 ELSE
512 cpassert(PRESENT(k_bf_start_from_atom))
513 my_k_offset = k_bf_start_from_atom(at_k) - 1
514 END IF
515 IF (PRESENT(atom_i)) THEN
516 my_i_offset = 0
517 ELSE
518 cpassert(PRESENT(i_bf_start_from_atom))
519 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 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 CALL gw_3c_ws_release(ws)
532 CALL gw_3c_ctx_release(ctx)
533
534 CALL timestop(handle)
535
536 END SUBROUTINE build_3c_integral_block
537
538END MODULE gw_integrals
539
imax
static int imax(int x, int y)
Returns the larger of two given integers (missing from the C standard)
Definition dbm_internal.h:14
cell_types::pbc
Definition cell_types.F:103
ai_contraction_sphi
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
Definition ai_contraction_sphi.F:15
ai_contraction_sphi::abc_contract_xsmm
subroutine, public abc_contract_xsmm(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer, prefac, pstfac)
3-center contraction routine from primitive cartesian Gaussians to spherical Gaussian functions using...
Definition ai_contraction_sphi.F:250
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:210
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::init_md_ftable
subroutine, public init_md_ftable(nmax)
Initialize a table of F_n(t) values in the range 0 <= t <= 12 with a stepsize of 0....
Definition gamma.F:540
gw_integrals
Utility method to build 3-center integrals for small cell GW.
Definition gw_integrals.F:11
gw_integrals::build_3c_integral_block
subroutine, public build_3c_integral_block(int_3c, qs_env, potential_parameter, basis_j, basis_k, basis_i, cell_j, cell_k, cell_i, atom_j, atom_k, atom_i, j_bf_start_from_atom, k_bf_start_from_atom, i_bf_start_from_atom)
...
Definition gw_integrals.F:461
gw_integrals::build_3c_integral_block_ctx
subroutine, public build_3c_integral_block_ctx(int_3c, ctx, ws, atom_j, atom_k, atom_i, cell_j, cell_k, cell_i, j_offset, k_offset, i_offset, screened)
Computes the 3c integral block (mu(atom_j) nu(atom_k) | P(atom_i)) for ONE atom triple,...
Definition gw_integrals.F:305
gw_integrals::gw_3c_ws_create
subroutine, public gw_3c_ws_create(ws, ctx)
Creates a per-thread 3c workspace: libint object + LIBXSMM contraction buffers.
Definition gw_integrals.F:259
gw_integrals::gw_3c_ws_release
subroutine, public gw_3c_ws_release(ws)
Releases a per-thread 3c workspace.
Definition gw_integrals.F:275
gw_integrals::gw_3c_ctx_release
subroutine, public gw_3c_ctx_release(ctx)
Releases the shared 3c-integral context.
Definition gw_integrals.F:232
gw_integrals::gw_3c_ctx_create
subroutine, public gw_3c_ctx_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i)
Builds the shared 3c-integral context: screening radii, basis maxima, contracted sphi tables,...
Definition gw_integrals.F:107
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:833
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:833
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:833
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:833
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_2c_3c::eri_3center
subroutine, public eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, int_abc_ext)
Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian gaussian orbit...
Definition libint_2c_3c.F:115
libint_2c_3c::cutoff_screen_factor
real(kind=dp), parameter, public cutoff_screen_factor
Definition libint_2c_3c.F:49
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
libint_wrapper::cp_libint_init_3eri
subroutine, public cp_libint_init_3eri(lib, max_am)
Definition libint_wrapper.F:1139
libint_wrapper::cp_libint_cleanup_3eri
subroutine, public cp_libint_cleanup_3eri(lib)
Definition libint_wrapper.F:1207
libint_wrapper::cp_libint_set_contrdepth
subroutine, public cp_libint_set_contrdepth(lib, contrdepth)
Definition libint_wrapper.F:1109
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:57
t_c_g0::init
subroutine, public init(nder, iunit, mepos, group)
...
Definition t_c_g0.F:1357
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition t_c_g0.F:1464
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_array_utils::cp_2d_r_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:109
gw_integrals::gw_3c_ctx_type
Shared read-only context for repeated 3-center integral block builds: screening parameters,...
Definition gw_integrals.F:65
gw_integrals::gw_3c_ws_type
Per-thread workspace for 3-center integral block builds: libint object + contraction buffers....
Definition gw_integrals.F:89
libint_2c_3c::libint_potential_type
Definition libint_2c_3c.F:67
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:71
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220