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 Interface to the LibGint-Library.
10 : !> \par History
11 : !> 10.2024 Created
12 : !> \author Marcello Puligheddu
13 : ! **************************************************************************************************
14 : MODULE libGint_wrapper
15 :
16 : USE kinds, ONLY: dp
17 : #if(__LIBGINT)
18 : USE input_constants, ONLY: do_potential_coulomb, do_potential_truncated
19 : USE libGint, ONLY: libgint_init, libgint_set_Potential_Truncated, libgint_set_hf_fac, libgint_set_max_mem, &
20 : libgint_set_P, libgint_set_P_polarized, libgint_set_K, libgint_set_K_polarized, &
21 : libgint_get_K, libgint_get_K_polarized, libgint_set_Atom, libgint_set_Atom_L, &
22 : libgint_set_cell, libgint_set_neighs, &
23 : libgint_add_prm, libgint_add_shell, libgint_add_cell, libgint_add_qrt, &
24 : libgint_add_qrtt, libgint_add_set
25 : USE t_c_g0, ONLY: C0
26 : #endif
27 : USE hfx_types, ONLY: hfx_type, hfx_memory_type, hfx_potential_type, &
28 : hfx_screen_coeff_type, hfx_cell_type, hfx_basis_type
29 :
30 : USE cell_types, ONLY: cell_type
31 : USE hfx_pair_list_methods, ONLY: build_pair_list_pbc_pgf, bra_t, allocate_bra
32 : USE particle_types, ONLY: particle_type
33 :
34 : USE iso_C_binding, ONLY: c_ptr
35 :
36 : #include "./base/base_uses.f90"
37 : IMPLICIT NONE
38 : PRIVATE
39 : INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: first_set_of_atom
40 : TYPE(bra_t), TARGET, SAVE :: bra, ket
41 : LOGICAL, SAVE :: first_call = .TRUE.
42 : TYPE(c_ptr), SAVE :: libGint_handle
43 : !$OMP THREADPRIVATE( first_set_of_atom, first_call )
44 : !$OMP THREADPRIVATE( bra,ket )
45 : !$OMP THREADPRIVATE( libGint_handle )
46 :
47 : PUBLIC :: cp_libGint_init, libGint_update_env, libGint_set_density, libGint_coulomb4, &
48 : libGint_update_fock_matrix, libGint_get_fock_matrix
49 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libGint_wrapper'
50 :
51 : ! Comunicates the current density to the libGint engine.
52 : ! Sets the Fock matrix on the device to zero
53 : INTERFACE libGint_set_density
54 : MODULE PROCEDURE libGint_set_density_A
55 : MODULE PROCEDURE libGint_set_density_AB
56 : END INTERFACE
57 :
58 : INTERFACE libGint_get_fock_matrix
59 : MODULE PROCEDURE libGint_get_fock_matrix_A
60 : MODULE PROCEDURE libGint_get_fock_matrix_AB
61 : END INTERFACE
62 :
63 : CONTAINS
64 : ! **************************************************************************************************
65 : !> \brief Sets libGint internal enviroment, must be called at least once before libGint can be used
66 : !> \param[in] actual_x_data pointer to hfx_data
67 : ! **************************************************************************************************
68 0 : SUBROUTINE cp_libGint_init(actual_x_data)
69 :
70 : TYPE(hfx_type), POINTER, INTENT(in) :: actual_x_data
71 : #if(__LIBGINT)
72 : ! Init the offload library, creates an handle unique to this omp thread
73 : CALL libgint_init(libGint_handle)
74 : ! Comunicate the chosen potential and its parameters to libGint
75 : IF (actual_x_data%potential_parameter%potential_type == do_potential_truncated) THEN
76 : ! truncated coulomb needs the C0 coefficients. We do not read or compute them,
77 : ! they must be already saved in C0 from the t[runcated]_c[oulomb]_g0 module
78 : CALL libgint_set_Potential_Truncated(libGint_handle, &
79 : actual_x_data%potential_parameter%cutoff_radius, &
80 : C0(:, :))
81 : ELSE
82 : CPABORT("The selected interaction potential type is not available with libGint")
83 : END IF
84 : first_call = .FALSE.
85 : #else
86 : MARK_USED(actual_x_data)
87 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
88 : #endif
89 0 : END SUBROUTINE cp_libGint_init
90 : ! **************************************************************************************************
91 : !> \brief Initialize and update the libGint computational environment. Must be called at least once
92 : !> after geo_change and before libGint can be used.
93 : !>
94 : !> This routine sets up data required by the libGint integral engine, including
95 : !> Hartree–Fock scaling factors, memory limits, periodic cell information, and
96 : !> per-atom/basis-set data. It also allocates CPU-side buffers (`bra` and `ket`)
97 : !> used for storing screened Gaussian primitive pairs.
98 : !>
99 : !> \param[in] fac Fraction of exact exchange
100 : !> \param[in] memory_parameter Pointer to memory configuration
101 : !> \param[in] do_periodic Logical flag: whether to consider pbc
102 : !> \param[in] cell primitive simulation cell
103 : !> \param[in] actual_x_data HF exchange data
104 : !> \param[in] nneighbors Lattice points
105 : !> \param[in] max_pgf Maximum number of primitive Gaussians per shel
106 : !> \param[in] natom Number of atoms in the system
107 : !> \param[in] kind_of Array mapping atom indices to atomic kinds
108 : !> \param[in] particle_set particle set, we extract the positions
109 : !> \param[in] basis_parameter gaussian basis set parameters
110 : !>
111 : ! **************************************************************************************************
112 0 : SUBROUTINE libGint_update_env(fac, memory_parameter, do_periodic, cell, actual_x_data, nneighbors, max_pgf, &
113 : natom, kind_of, particle_set, basis_parameter)
114 :
115 : REAL(dp) :: fac
116 : TYPE(hfx_memory_type), POINTER :: memory_parameter
117 : LOGICAL :: do_periodic
118 : TYPE(cell_type), POINTER :: cell
119 : TYPE(hfx_type), POINTER :: actual_x_data
120 : INTEGER :: nneighbors, max_pgf, natom
121 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
122 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
123 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
124 : #if(__LIBGINT)
125 : LOGICAL(1) :: do_pbc
126 : REAL(dp), DIMENSION(:, :), ALLOCATABLE :: cell_r
127 : REAL(dp), DIMENSION(:), ALLOCATABLE :: flat_gcc
128 : INTEGER, DIMENSION(:), POINTER :: la_min, la_max, npgfa, nsgfa
129 : INTEGER, DIMENSION(:, :), POINTER :: nsgfl_a
130 : INTEGER :: i, l, iset, jset, iatom, ikind, nseta, inla, nla
131 : REAL(dp) :: ra(3)
132 : REAL(dp), DIMENSION(:, :), POINTER :: zeta
133 : REAL(dp), DIMENSION(:, :, :), POINTER :: gcc
134 : TYPE(bra_t), POINTER :: bra_p, ket_p
135 :
136 : ! Set the multiplicative factor fac (the fraction of exact exchange times the spin factor) for libGint
137 : CALL libGint_set_hf_fac(libGint_handle, fac)
138 : ! Comunicate max gpu mem per mpi
139 : CALL libGint_set_max_mem(libGint_handle, memory_parameter%max_memory)
140 :
141 : ! Info about periodic cells and neighbouring cells
142 : do_pbc = do_periodic
143 : ALLOCATE (cell_r(3, nneighbors))
144 : DO i = 1, nneighbors
145 : cell_r(:, i) = actual_x_data%neighbor_cells(i)%cell_r(:)
146 : END DO
147 : CALL libgint_set_cell(libGint_handle, do_pbc, cell%hmat, cell%h_inv)
148 : CALL libgint_set_neighs(libGint_handle, cell_r, nneighbors)
149 :
150 : ! CPU side temporary arrays for info about the <AB(g) and CD(n)> pairs
151 : bra_p => bra
152 : ket_p => ket
153 : CALL allocate_bra(bra_p, max_pgf, nneighbors)
154 : CALL allocate_bra(ket_p, max_pgf, nneighbors)
155 :
156 : ! Comunicate atomset and atom info to the engine
157 : jset = 1
158 : IF (ALLOCATED(first_set_of_atom)) DEALLOCATE (first_set_of_atom)
159 : ALLOCATE (first_set_of_atom(natom))
160 : DO iatom = 1, natom
161 : ikind = kind_of(iatom)
162 : ra = particle_set(iatom)%r(:)
163 : nseta = basis_parameter(ikind)%nset
164 : npgfa => basis_parameter(ikind)%npgf
165 : la_min => basis_parameter(ikind)%lmin
166 : la_max => basis_parameter(ikind)%lmax
167 : zeta => basis_parameter(ikind)%zet
168 : nsgfa => basis_parameter(ikind)%nsgf
169 : nsgfl_a => basis_parameter(ikind)%nsgfl
170 : gcc => basis_parameter(ikind)%gcc
171 : first_set_of_atom(iatom) = jset
172 : DO iset = 1, nseta
173 : CALL libgint_set_Atom(libGint_handle, jset - 1, ra, zeta(:, iset), npgfa(iset))
174 : inla = 1
175 : DO l = la_min(iset), la_max(iset)
176 : nla = nsgfl_a(l, iset)
177 : IF (ALLOCATED(flat_gcc)) DEALLOCATE (flat_gcc)
178 : ALLOCATE (flat_gcc(npgfa(iset)*nla))
179 : flat_gcc = PACK(gcc(1:npgfa(iset), inla:inla + nla - 1, iset), .TRUE.)
180 : CALL libgint_set_Atom_L(libGint_handle, jset - 1, l, nla, flat_gcc)
181 : inla = inla + nla
182 : END DO
183 : jset = jset + 1
184 : END DO
185 : END DO
186 : #else
187 : MARK_USED(fac)
188 : MARK_USED(memory_parameter)
189 : MARK_USED(do_periodic)
190 : MARK_USED(cell)
191 : MARK_USED(actual_x_data)
192 : MARK_USED(nneighbors)
193 : MARK_USED(max_pgf)
194 : MARK_USED(natom)
195 : MARK_USED(kind_of)
196 : MARK_USED(particle_set)
197 : MARK_USED(basis_parameter)
198 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
199 : #endif
200 0 : END SUBROUTINE libGint_update_env
201 :
202 : ! **************************************************************************************************
203 : !> \brief communicates the density to libGint, no spin case
204 : !>
205 : !> \param[in] full_density_alpha full density matrix array
206 : ! **************************************************************************************************
207 0 : SUBROUTINE libGint_set_density_A(full_density_alpha)
208 :
209 : REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha
210 : #if(__LIBGINT)
211 : CALL libgint_set_P(libGint_handle, full_density_alpha(:, 1))
212 : #else
213 : MARK_USED(full_density_alpha)
214 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
215 : #endif
216 0 : END SUBROUTINE libGint_set_density_A
217 :
218 : ! **************************************************************************************************
219 : !> \brief communicates the density to libGint, spin case
220 : !>
221 : !> \param[in] full_density_alpha full density matrix array on the alpha channel
222 : !> \param[in] full_density_beta full density matrix array on the beta channel
223 : ! **************************************************************************************************
224 0 : SUBROUTINE libGint_set_density_AB(full_density_alpha, full_density_beta)
225 :
226 : REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha, full_density_beta
227 : #if(__LIBGINT)
228 : CALL libgint_set_P_polarized(libGint_handle, full_density_alpha(:, 1), full_density_beta(:, 1))
229 : #else
230 : MARK_USED(full_density_alpha)
231 : MARK_USED(full_density_beta)
232 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
233 : #endif
234 0 : END SUBROUTINE libGint_set_density_AB
235 :
236 : ! **************************************************************************************************
237 : !> \brief requests the Fock matrix ( hf_fraction * D @@ I ) from libGint. When this function
238 : !> returns, full_ks_alpha_from_gpu will contain the fock matrix for this MPI rank.
239 : !>
240 : !> \param[out] full_ks_alpha_from_gpu location, already allocated, for the Fock matrix
241 : !> \note assumes full_ks_alpha_from_gpu is already allocated with enough memory
242 : ! **************************************************************************************************
243 0 : SUBROUTINE libGint_get_fock_matrix_A(full_ks_alpha_from_gpu)
244 :
245 : REAL(dp), DIMENSION(:, :), POINTER :: full_ks_alpha_from_gpu
246 : #if(__LIBGINT)
247 : CALL libGint_get_K(libGint_handle, full_ks_alpha_from_gpu(:, 1))
248 : #else
249 : MARK_USED(full_ks_alpha_from_gpu)
250 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
251 : #endif
252 0 : END SUBROUTINE libGint_get_fock_matrix_A
253 :
254 : ! **************************************************************************************************
255 : !> \brief requests the Fock matrices from libGint for both channels. When this function
256 : !> returns, full_ks_alpha_from_gpu and beta will contain the fock matrix for this MPI rank.
257 : !>
258 : !> \param[out] full_ks_alpha_from_gpu location, already allocated, for the Fock matrix alpha
259 : !> \param[out] full_ks_beta_from_gpu location, already allocated, for the Fock matrix beta
260 : !> \note assumes full_ks_alpha_from_gpu and full_ks_beta_from_gpu are already allocated
261 : ! **************************************************************************************************
262 0 : SUBROUTINE libGint_get_fock_matrix_AB(full_ks_alpha_from_gpu, full_ks_beta_from_gpu)
263 :
264 : REAL(dp), DIMENSION(:, :), POINTER :: full_ks_alpha_from_gpu, full_ks_beta_from_gpu
265 : #if(__LIBGINT)
266 : CALL libgint_get_K_polarized(libGint_handle, full_ks_alpha_from_gpu(:, 1), full_ks_beta_from_gpu(:, 1))
267 : #else
268 : MARK_USED(full_ks_alpha_from_gpu)
269 : MARK_USED(full_ks_beta_from_gpu)
270 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
271 : #endif
272 0 : END SUBROUTINE libGint_get_fock_matrix_AB
273 :
274 : ! **************************************************************************************************
275 : !> \brief Assign two-electron integrals of a quartet/shell to libGint
276 : !>
277 : !> First we build a list of ab primitives and cd primitives using build_pair_list_pbc_pgf
278 : !> These lists will contain (before screening) n_cell * npgf * npgf.
279 : !> We loop over both and send ab(g) | cd(h) to libgint, which will loop internally over n
280 : !>
281 : !> \param[in] iatom index of atom I / A
282 : !> \param[in] jatom index of atom J / B
283 : !> \param[in] katom index of atom K / C
284 : !> \param[in] latom index of atom L / D
285 : !> \param[in] iset index of the set for atom iatom we are adding
286 : !> \param[in] jset index of the set for atom jatom we are adding
287 : !> \param[in] kset index of the set for atom katom we are adding
288 : !> \param[in] lset index of the set for atom latom we are adding
289 : !> \param[in] ra position of atom iatom
290 : !> \param[in] rb position of atom jatom
291 : !> \param[in] rc position of atom katom
292 : !> \param[in] rd position of atom latom
293 : !> \param[in] npgfa number of gaussians in set iset of atom iatom
294 : !> \param[in] npgfb number of gaussians in set jset of atom jatom
295 : !> \param[in] npgfc number of gaussians in set kset of atom katom
296 : !> \param[in] npgfd number of gaussians in set lset of atom latom
297 : !> \param[in] potential_parameter information about the potential
298 : !> \param[in] screen1 screening information for AB pair
299 : !> \param[in] screen2 screening information for CD pair
300 : !> \param[in] log10_pmax density screening factor
301 : !> \param[in] log10_eps_schwarz screening tolerance
302 : !> \param[in] pgf1 screening information for each AB pair primitive
303 : !> \param[in] pgf2 screening information for each CD pair primitive
304 : !> \param[in] neighbor_cells array with lattice vectors
305 : !> \param[in] cell information about the simulation box
306 : !> \param[in] do_periodic flag for pbc
307 : !> \param[out] screened whether the whole quartet was screened out
308 : ! **************************************************************************************************
309 :
310 0 : SUBROUTINE libGint_coulomb4(iatom, jatom, katom, latom, iset, jset, kset, lset, &
311 : ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, &
312 : potential_parameter, &
313 : screen1, screen2, log10_pmax, log10_eps_schwarz, &
314 : pgf1, pgf2, &
315 : neighbor_cells, cell, do_periodic, screened)
316 :
317 : INTEGER, INTENT(in) :: iatom, jatom, katom, latom, iset, jset, kset, lset
318 : REAL(dp), INTENT(IN) :: ra(3), rb(3), rc(3), rd(3)
319 : INTEGER, INTENT(IN) :: npgfa, npgfb, npgfc, npgfd
320 : TYPE(hfx_potential_type) :: potential_parameter
321 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2)
322 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
323 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
324 : POINTER :: pgf1, pgf2
325 : TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
326 : TYPE(cell_type), POINTER :: cell
327 : LOGICAL, INTENT(IN) :: do_periodic
328 : LOGICAL, INTENT(out) :: screened
329 :
330 : #if(__LIBGINT)
331 : TYPE(bra_t), POINTER :: bra_p, ket_p
332 : LOGICAL :: cell_was_screened
333 : INTEGER :: idx_n1, idx_n2, n1, n2, idx_ij, idx_kl, o_ij, n_ij, o_kl, n_kl !, n3
334 : INTEGER :: ipgf, jpgf, kpgf, lpgf, iatom_set, jatom_set, katom_set, latom_set
335 : REAL(dp) :: pgf_max_1, pgf_max_2 ! , R1, R2, rpq2
336 : INTEGER :: nelements_ij, nelements_kl
337 :
338 : cell = cell
339 : potential_parameter = potential_parameter
340 : bra_p => bra
341 : ket_p => ket
342 :
343 : screened = .TRUE.
344 : iatom_set = first_set_of_atom(iatom) + iset - 2
345 : jatom_set = first_set_of_atom(jatom) + jset - 2
346 : katom_set = first_set_of_atom(katom) + kset - 2
347 : latom_set = first_set_of_atom(latom) + lset - 2
348 :
349 : CALL build_pair_list_pbc_pgf(npgfa, npgfb, bra_p, screen1, screen2, &
350 : pgf1, log10_pmax, log10_eps_schwarz, ra, rb, &
351 : nelements_ij, neighbor_cells, do_periodic)
352 :
353 : CALL build_pair_list_pbc_pgf(npgfc, npgfd, ket_p, screen2, screen1, &
354 : pgf2, log10_pmax, log10_eps_schwarz, rc, rd, &
355 : nelements_kl, neighbor_cells, do_periodic)
356 :
357 : ! Note: we use 3 numbers n1 n2 and n3 as indices for the lattice traslantion vectors
358 : ! n1 for the AB pair, n2 for the CD pair and n3 for the PQ pair
359 : ! so that e.g. B = B0 + ua(n1) means B.x = B0.x + ua(n1).x (and same for y and z)
360 : ! the ua, saved in neighbor_cells(:)%cell_r(:), are already computed
361 : ! lattice translation vectors T = i a1 + j a2 + k a3
362 : ! where a1,a2 and a3 are lattice vectors and i,j and k integers.
363 : ! So, B.x = B0.x + n1.i * a1.x + n1.j * a2.x + n1.k * a3.x (and same for y and z)
364 : !
365 :
366 : DO idx_n1 = 1, bra%cell_cnt
367 :
368 : n1 = bra%cell_idx(1, idx_n1)
369 : n_ij = bra%cell_idx(2, idx_n1)
370 : o_ij = bra%cell_idx(3, idx_n1)
371 :
372 : DO idx_n2 = 1, ket%cell_cnt
373 :
374 : n2 = ket%cell_idx(1, idx_n2)
375 : n_kl = ket%cell_idx(2, idx_n2)
376 : o_kl = ket%cell_idx(3, idx_n2)
377 :
378 : cell_was_screened = .TRUE.
379 : DO idx_ij = o_ij + 1, o_ij + n_ij
380 :
381 : ipgf = bra%pgf_idx(1, idx_ij)
382 : jpgf = bra%pgf_idx(2, idx_ij)
383 :
384 : pgf_max_1 = bra%pgf_scr(1, idx_ij)
385 :
386 : DO idx_kl = o_kl + 1, o_kl + n_kl
387 : kpgf = ket%pgf_idx(1, idx_kl)
388 : lpgf = ket%pgf_idx(2, idx_kl)
389 :
390 : pgf_max_2 = ket%pgf_scr(1, idx_kl)
391 :
392 : IF (pgf_max_1 + pgf_max_2 + log10_pmax < log10_eps_schwarz) CYCLE
393 :
394 : CALL libGint_add_prm(libGint_handle, ipgf - 1, jpgf - 1, kpgf - 1, lpgf - 1)
395 : cell_was_screened = .FALSE.
396 :
397 : END DO ! ket pgf
398 : END DO ! bra pgf
399 :
400 : IF (.NOT. cell_was_screened) THEN
401 : CALL libgint_add_shell(libGint_handle, iatom_set, jatom_set, katom_set, latom_set, n1 - 1, n2 - 1)
402 : screened = .FALSE.
403 : END IF
404 :
405 : END DO ! ket n2
406 : END DO ! bra n1
407 : #else
408 : MARK_USED(iatom)
409 : MARK_USED(jatom)
410 : MARK_USED(katom)
411 : MARK_USED(latom)
412 : MARK_USED(iset)
413 : MARK_USED(jset)
414 : MARK_USED(kset)
415 : MARK_USED(lset)
416 : MARK_USED(ra)
417 : MARK_USED(rb)
418 : MARK_USED(rc)
419 : MARK_USED(rd)
420 : MARK_USED(npgfa)
421 : MARK_USED(npgfb)
422 : MARK_USED(npgfc)
423 : MARK_USED(npgfd)
424 : MARK_USED(potential_parameter)
425 : MARK_USED(screen1)
426 : MARK_USED(screen2)
427 : MARK_USED(log10_pmax)
428 : MARK_USED(log10_eps_schwarz)
429 : MARK_USED(pgf1)
430 : MARK_USED(pgf2)
431 : MARK_USED(neighbor_cells)
432 : MARK_USED(cell)
433 : MARK_USED(do_periodic)
434 : MARK_USED(screened)
435 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
436 : #endif
437 0 : END SUBROUTINE libGint_coulomb4
438 :
439 : ! **************************************************************************************************
440 : !> \brief The previous coulomb_4 function assigned an integral between primitive gaussian.
441 : !> Now we assign the nsgfa(iset)*b*c*d gcc integrals form this set,
442 : !> along with the Pac Pad Pbc Pbd density to the Kbd Kbc Kad Kac Fock matrix
443 : !>
444 : !> \param[in] symm_fac simmetry factor from iatom=jatom and the like
445 : !> \param[in] iatom index of atom A
446 : !> \param[in] jatom index of atom B
447 : !> \param[in] katom index of atom C
448 : !> \param[in] latom index of atom D
449 : !> \param[in] iset index of set for atom A
450 : !> \param[in] jset index of set for atom B
451 : !> \param[in] kset index of set for atom C
452 : !> \param[in] lset index of set for atom D
453 : !> \param[in] atomic_offset_ac global offset for the pair of A and C atom in the density ( and Fock) matrix
454 : !> \param[in] atomic_offset_ad global offset for the pair of A and D atom in the density ( and Fock) matrix
455 : !> \param[in] atomic_offset_bc global offset for the pair of B and C atom in the density ( and Fock) matrix
456 : !> \param[in] atomic_offset_bd global offset for the pair of B and D atom in the density ( and Fock) matrix
457 : !> \param[in] offset_ac_set matrix of sub_offset for sets in atomic_offset_ac
458 : !> \param[in] offset_ad_set matrix of sub_offset for sets in atomic_offset_ad
459 : !> \param[in] offset_bc_set matrix of sub_offset for sets in atomic_offset_bc
460 : !> \param[in] offset_bd_set matrix of sub_offset for sets in atomic_offset_bd
461 : !> \param[in] nsgfa total number of (spherical, contracted) integrals for iset.
462 : !> Used as leading dimension of the AC, AD sublock, depending on transposition
463 : !> \param[in] nsgfb total number of (spherical, contracted) integrals for jset.
464 : !> Used as leading dimension of the BC, BD sublock, depending on transposition
465 : !> \param[in] nsgfc total number of (spherical, contracted) integrals for kset.
466 : !> Used as leading dimension of the AC, BC sublock, depending on transposition
467 : !> \param[in] nsgfd total number of (spherical, contracted) integrals for lset.
468 : !> Used as leading dimension of the AD, BD sublock, depending on transposition
469 : !> \param[in] la_min minumum angular moment in set iset of atom A
470 : !> \param[in] la_max maximum angular moment in set iset of atom A
471 : !> \param[in] lb_min minumum angular moment in set iset of atom B
472 : !> \param[in] lb_max maximum angular moment in set iset of atom B
473 : !> \param[in] lc_min minumum angular moment in set iset of atom C
474 : !> \param[in] lc_max maximum angular moment in set iset of atom C
475 : !> \param[in] ld_min minumum angular moment in set iset of atom D
476 : !> \param[in] ld_max maximum angular moment in set iset of atom D
477 : !> \param[in] nsgfl_a matrix with the number of linear combinations of primitive gaussians
478 : !> for each angular moment for each set in atom A. We only read iset
479 : !> \param[in] nsgfl_b matrix with the number of linear combinations of primitive gaussians
480 : !> for each angular moment for each set in atom B. We only read jset
481 : !> \param[in] nsgfl_c matrix with the number of linear combinations of primitive gaussians
482 : !> for each angular moment for each set in atom C. We only read kset
483 : !> \param[in] nsgfl_d matrix with the number of linear combinations of primitive gaussians
484 : !> for each angular moment for each set in atom D. We only read lset
485 : !> \note
486 : !> The atomic_offset_xy, offset_xy_set matrices provide set offsets which are combined with per-L and
487 : !> per linear combination offsets to produce the final indices into the dense (but block-sparse)
488 : !> density and Fock matrices
489 : !> - The routine assumes Fortran column major order and contiguous storage for the per set
490 : !> density subblocks as described in the code comments.
491 : !> - The code computes transposition flags to only use the lower part of P and K
492 : !>
493 : !>
494 : ! **************************************************************************************************
495 :
496 0 : SUBROUTINE libGint_update_fock_matrix( &
497 : symm_fac, &
498 : iatom, jatom, katom, latom, &
499 : iset, jset, kset, lset, &
500 : atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd, &
501 : offset_ac_set, offset_ad_set, offset_bc_set, offset_bd_set, &
502 : nsgfa, nsgfb, nsgfc, nsgfd, &
503 : la_min, la_max, lb_min, lb_max, &
504 : lc_min, lc_max, ld_min, ld_max, &
505 : nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d)
506 :
507 : REAL(dp) :: symm_fac
508 : INTEGER :: iatom, jatom, katom, latom
509 : INTEGER :: iset, jset, kset, lset
510 : INTEGER :: atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd
511 : INTEGER, DIMENSION(:, :), POINTER :: offset_ac_set, offset_ad_set
512 : INTEGER, DIMENSION(:, :), POINTER :: offset_bc_set, offset_bd_set
513 : INTEGER :: nsgfa, nsgfb, nsgfc, nsgfd
514 : INTEGER :: la_min, la_max, lb_min, lb_max
515 : INTEGER :: lc_min, lc_max, ld_min, ld_max
516 : INTEGER, DIMENSION(:, :), POINTER :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d
517 :
518 : #if(__LIBGINT)
519 : !! (Hyp)
520 : ! Let a be a set composed of 2 s and 1 p function.
521 : ! Let c be a set composed of 1 s and 2 p function.
522 : ! (1) The density matrix for the ac pair is a 5 x 7 matrix organized as
523 : !
524 : ! / -------------------------------------------------------------------------------------------------------------------\
525 : ! | a_s1_0@c_s1_0 || a_s1_0@c_p1_0 | a_s1_0@c_p1_1 | a_s1_0@c_p1_2 || a_s1_0@c_p2_0 | a_s1_0@c_p2_1 | a_s1_0@c_p2_2 |
526 : ! | a_s2_0@c_s1_0 || a_s2_0@c_p1_0 | a_s2_0@c_p1_1 | a_s2_0@c_p1_2 || a_s2_0@c_p2_0 | a_s2_0@c_p2_1 | a_s2_0@c_p2_2 |
527 : ! | a_p1_0@c_s1_0 || a_p1_0@c_p1_0 | a_p1_0@c_p1_1 | a_p1_0@c_p1_2 || a_p1_0@c_p2_0 | a_p1_0@c_p2_1 | a_p1_0@c_p2_2 |
528 : ! | a_p1_1@c_s1_0 || a_p1_1@c_p1_0 | a_p1_1@c_p1_1 | a_p1_1@c_p1_2 || a_p1_1@c_p2_0 | a_p1_1@c_p2_1 | a_p1_1@c_p2_2 |
529 : ! | a_p1_2@c_s1_0 || a_p1_2@c_p1_0 | a_p1_2@c_p1_1 | a_p1_2@c_p1_2 || a_p1_2@c_p2_0 | a_p1_2@c_p2_1 | a_p1_2@c_p2_2 |
530 : ! \ -------------------------------------------------------------------------------------------------------------------/
531 : !
532 : ! where A_LX_Y means the (Y+1) component of the Xth linear combination of the L angular moment for atom A
533 : !
534 : ! (2) This matrix is dense, rectangular and contigous in memory, in fortran column major order.
535 : ! (3) The big matrix with all pairs is block sparse triangular, only the lower part is valid.
536 :
537 : LOGICAL(1) :: Tac, Tad, Tbc, Tbd
538 : INTEGER :: offset_ac_L_set, offset_ad_L_set, offset_bc_L_set, offset_bd_L_set
539 : INTEGER :: s_offset_a, s_offset_b, s_offset_c, s_offset_d
540 : INTEGER :: s_offset_a_l, s_offset_b_l, s_offset_c_l, s_offset_d_l
541 : INTEGER :: s_offset_ac, s_offset_ad, s_offset_bc, s_offset_bd
542 : INTEGER :: ld_ac_set, ld_ad_set, ld_bc_set, ld_bd_set
543 : INTEGER :: la, lb, lc, ld, nla, nlb, nlc, nld, inla, inlb, inlc, inld
544 : ! TODO rewrite as update_fock_matrix_gpu(libGint_handle, iatomset,jatomset,katomset,latomset )
545 : ! AFTER TODO communicate (the pointer to) atomic_offset to libGint AND
546 : ! AFTER TODO communicate (the pointer to) set_offset to libGint AND
547 : ! AFTER TODO check if this idea makes sense in general for other codes
548 : !
549 : ! Note: this would not change the need to compute sub offsets and
550 : ! the 8 loops, it would just transfer them to libGint
551 : ! Except, if libGint can be sure every set has 1 l, it can collapse the l loops
552 : ! and/or, if libGint can be sure every l has 1 nl, it can collapse the n loops
553 : IF (jatom >= latom) THEN
554 : offset_bd_L_set = offset_bd_set(jset, lset) + atomic_offset_bd - 2
555 : ld_bd_set = nsgfb
556 : Tbd = .FALSE.
557 : ELSE
558 : offset_bd_L_set = offset_bd_set(lset, jset) + atomic_offset_bd - 2
559 : ld_bd_set = nsgfd
560 : Tbd = .TRUE.
561 : END IF
562 : IF (jatom >= katom) THEN
563 : offset_bc_L_set = offset_bc_set(jset, kset) + atomic_offset_bc - 2
564 : ld_bc_set = nsgfb
565 : Tbc = .FALSE.
566 : ELSE
567 : offset_bc_L_set = offset_bc_set(kset, jset) + atomic_offset_bc - 2
568 : ld_bc_set = nsgfc
569 : Tbc = .TRUE.
570 : END IF
571 :
572 : IF (iatom >= latom) THEN
573 : offset_ad_L_set = offset_ad_set(iset, lset) + atomic_offset_ad - 2
574 : ld_ad_set = nsgfa
575 : Tad = .FALSE.
576 : ELSE
577 : offset_ad_L_set = offset_ad_set(lset, iset) + atomic_offset_ad - 2
578 : ld_ad_set = nsgfd
579 : Tad = .TRUE.
580 : END IF
581 :
582 : IF (iatom >= katom) THEN
583 : offset_ac_L_set = offset_ac_set(iset, kset) + atomic_offset_ac - 2
584 : ld_ac_set = nsgfa
585 : Tac = .FALSE.
586 : ELSE
587 : offset_ac_L_set = offset_ac_set(kset, iset) + atomic_offset_ac - 2
588 : ld_ac_set = nsgfc
589 : Tac = .TRUE.
590 : END IF
591 :
592 : s_offset_a_l = 0
593 : DO la = la_min, la_max
594 : nla = nsgfl_a(la, iset)
595 : s_offset_b_l = 0
596 : DO lb = lb_min, lb_max
597 : nlb = nsgfl_b(lb, jset)
598 : s_offset_c_l = 0
599 : DO lc = lc_min, lc_max
600 : nlc = nsgfl_c(lc, kset)
601 : s_offset_d_l = 0
602 : ld_loop: DO ld = ld_min, ld_max
603 : nld = nsgfl_d(ld, lset)
604 : CALL libgint_add_qrt(libGint_handle, la, lb, lc, ld, nla, nlb, nlc, nld)
605 : DO inla = 1, nla
606 : s_offset_a = s_offset_a_l + (inla - 1)*(2*la + 1)
607 : DO inlb = 1, nlb
608 : s_offset_b = s_offset_b_l + (inlb - 1)*(2*lb + 1)
609 : DO inlc = 1, nlc
610 : s_offset_c = s_offset_c_l + (inlc - 1)*(2*lc + 1)
611 : DO inld = 1, nld
612 : s_offset_d = s_offset_d_l + (inld - 1)*(2*ld + 1)
613 : IF (.NOT. Tac) THEN
614 : s_offset_ac = offset_ac_L_set + s_offset_c*ld_ac_set + s_offset_a
615 : ELSE
616 : s_offset_ac = offset_ac_L_set + s_offset_a*ld_ac_set + s_offset_c
617 : END IF
618 :
619 : IF (.NOT. Tad) THEN
620 : s_offset_ad = offset_ad_L_set + s_offset_d*ld_ad_set + s_offset_a
621 : ELSE
622 : s_offset_ad = offset_ad_L_set + s_offset_a*ld_ad_set + s_offset_d
623 : END IF
624 :
625 : IF (.NOT. Tbc) THEN
626 : s_offset_bc = offset_bc_L_set + s_offset_c*ld_bc_set + s_offset_b
627 : ELSE
628 : s_offset_bc = offset_bc_L_set + s_offset_b*ld_bc_set + s_offset_c
629 : END IF
630 :
631 : IF (.NOT. Tbd) THEN
632 : s_offset_bd = offset_bd_L_set + s_offset_d*ld_bd_set + s_offset_b
633 : ELSE
634 : s_offset_bd = offset_bd_L_set + s_offset_b*ld_bd_set + s_offset_d
635 : END IF
636 :
637 : CALL libgint_add_qrtt(libGint_handle, symm_fac, &
638 : la, lb, lc, ld, inla - 1, inlb - 1, inlc - 1, inld - 1, &
639 : ld_ac_set, ld_ad_set, ld_bc_set, ld_bd_set, &
640 : s_offset_ac, s_offset_ad, s_offset_bc, s_offset_bd, &
641 : Tac, Tad, Tbc, Tbd)
642 :
643 : END DO
644 : END DO
645 : END DO
646 : END DO
647 : s_offset_d_l = s_offset_d_l + nld*(2*ld + 1)
648 : END DO ld_loop
649 : s_offset_c_l = s_offset_c_l + nlc*(2*lc + 1)
650 : END DO
651 : s_offset_b_l = s_offset_b_l + nlb*(2*lb + 1)
652 : END DO
653 : s_offset_a_l = s_offset_a_l + nla*(2*la + 1)
654 : END DO
655 :
656 : CALL libgint_add_set(libGint_handle)
657 :
658 : #else
659 : MARK_USED(symm_fac)
660 : MARK_USED(iatom)
661 : MARK_USED(jatom)
662 : MARK_USED(katom)
663 : MARK_USED(latom)
664 : MARK_USED(iset)
665 : MARK_USED(jset)
666 : MARK_USED(kset)
667 : MARK_USED(lset)
668 : MARK_USED(atomic_offset_ac)
669 : MARK_USED(atomic_offset_ad)
670 : MARK_USED(atomic_offset_bc)
671 : MARK_USED(atomic_offset_bd)
672 : MARK_USED(offset_ac_set)
673 : MARK_USED(offset_ad_set)
674 : MARK_USED(offset_bc_set)
675 : MARK_USED(offset_bd_set)
676 : MARK_USED(nsgfa)
677 : MARK_USED(nsgfb)
678 : MARK_USED(nsgfc)
679 : MARK_USED(nsgfd)
680 : MARK_USED(la_min)
681 : MARK_USED(la_max)
682 : MARK_USED(lb_min)
683 : MARK_USED(lb_max)
684 : MARK_USED(lc_min)
685 : MARK_USED(lc_max)
686 : MARK_USED(ld_min)
687 : MARK_USED(ld_max)
688 : MARK_USED(nsgfl_a)
689 : MARK_USED(nsgfl_b)
690 : MARK_USED(nsgfl_c)
691 : MARK_USED(nsgfl_d)
692 0 : CPABORT("This CP2K executable has not been linked against the required library libGint.")
693 : #endif
694 0 : END SUBROUTINE libGint_update_fock_matrix
695 :
696 : END MODULE libGint_wrapper
|
