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