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 Contraction of integrals over primitive Cartesian Gaussians based on the contraction
10 : !> matrix sphi which is part of the gto_basis_set_type
11 : !> \par History
12 : !> -added abc_contract_xsmm routine, A. Bussy (04.2020)
13 : !> \author Dorothea Golze (05.2016)
14 : ! **************************************************************************************************
15 : MODULE ai_contraction_sphi
16 :
17 : USE kinds, ONLY: dp, int_8
18 : #if defined(__LIBXS)
19 : USE LIBXS_JIT, ONLY: libxs_gemm_config_t, libxs_gemm_dispatch, &
20 : libxs_gemm_call, LIBXS_DATATYPE_F64, C_LOC
21 : #endif
22 :
23 : #include "../base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 :
27 : PRIVATE
28 :
29 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_contraction_sphi'
30 :
31 : PUBLIC :: ab_contract, abc_contract, abcd_contract, abc_contract_xsmm
32 :
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief contract overlap integrals (a,b) and transfer to spherical Gaussians
37 : !> \param abint contracted, normalized integrals of spherical Gaussians
38 : !> \param sab uncontracted, unnormalized integrals of primitive Cartesian Gaussians
39 : !> \param sphi_a contraction matrix for center a
40 : !> \param sphi_b contraction matrix for center b
41 : !> \param ncoa number of cartesian orbitals on a
42 : !> \param ncob number of cartesian orbitals on b
43 : !> \param nsgfa number of spherical Gaussian functions on a
44 : !> \param nsgfb number of spherical Gaussian functions on b
45 : ! **************************************************************************************************
46 1169586 : SUBROUTINE ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
47 :
48 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: abint
49 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sab, sphi_a, sphi_b
50 : INTEGER, INTENT(IN) :: ncoa, ncob, nsgfa, nsgfb
51 :
52 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ab_contract'
53 :
54 : INTEGER :: handle
55 1169586 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cpp
56 :
57 1169586 : CALL timeset(routineN, handle)
58 :
59 1169586 : CPASSERT(ncob <= SIZE(sab, 2))
60 :
61 1169586 : IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (nsgfb*ncoa*(ncob + nsgfa))) THEN ! (sphi_a x sab) x sphi_b
62 3186868 : ALLOCATE (cpp(nsgfa, ncob))
63 : ! [nsgfa,ncoa] x [ncoa,ncob] -> [nsgfa,ncob]
64 24630064 : CALL dgemm("T", "N", nsgfa, ncob, ncoa, 1._dp, sphi_a, SIZE(sphi_a, 1), sab, SIZE(sab, 1), 0.0_dp, cpp, nsgfa)
65 : ! [nsgfa,ncob] x [ncob,nsgfb] -> [nsgfa,nsgfb]
66 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpp, nsgfa, sphi_b, SIZE(sphi_b, 1), 0.0_dp, &
67 54550439 : abint, SIZE(abint, 1))
68 : ELSE ! sphi_a x (sab x sphi_b)
69 1491476 : ALLOCATE (cpp(ncoa, nsgfb))
70 : ! [ncoa,ncob] x [ncob,nsgfb] -> [ncoa,nsgfb]
71 8303829 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1._dp, sab, SIZE(sab, 1), sphi_b, SIZE(sphi_b, 1), 0.0_dp, cpp, ncoa)
72 : ! [nsgfa,ncoa] x [ncoa,nsgfb] -> [nsgfa,nsgfb]
73 : CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1._dp, sphi_a, SIZE(sphi_a, 1), cpp, ncoa, 0.0_dp, &
74 13928925 : abint, SIZE(abint, 1))
75 : END IF
76 :
77 1169586 : DEALLOCATE (cpp)
78 :
79 1169586 : CALL timestop(handle)
80 :
81 1169586 : END SUBROUTINE ab_contract
82 :
83 : ! **************************************************************************************************
84 : !> \brief contract three-center overlap integrals (a,b,c) and transfer
85 : !> to spherical Gaussians
86 : !> \param abcint contracted, normalized integrals of spherical Gaussians
87 : !> \param sabc uncontracted, unnormalized integrals of primitive Cartesian Gaussians
88 : !> \param sphi_a contraction matrix for center a
89 : !> \param sphi_b contraction matrix for center b
90 : !> \param sphi_c contraction matrix for center c
91 : !> \param ncoa number of cartesian orbitals on a
92 : !> \param ncob number of cartesian orbitals on b
93 : !> \param ncoc number of cartesian orbitals on c
94 : !> \param nsgfa number of spherical Gaussian functions on a
95 : !> \param nsgfb number of spherical Gaussian functions on b
96 : !> \param nsgfc number of spherical Gaussian functions on c
97 : ! **************************************************************************************************
98 126004 : SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
99 : nsgfa, nsgfb, nsgfc)
100 :
101 : REAL(KIND=dp), DIMENSION(:, :, :) :: abcint, sabc
102 : REAL(KIND=dp), DIMENSION(:, :) :: sphi_a, sphi_b, sphi_c
103 : INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
104 :
105 : CHARACTER(LEN=*), PARAMETER :: routineN = 'abc_contract'
106 :
107 : INTEGER :: handle, i, m1, m2, m3, msphia, msphib, &
108 : msphic, mx
109 126004 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: tmp
110 :
111 126004 : CALL timeset(routineN, handle)
112 :
113 126004 : CPASSERT(SIZE(abcint, 1) == nsgfa)
114 126004 : CPASSERT(SIZE(abcint, 2) == nsgfb)
115 :
116 126004 : msphia = SIZE(sphi_a, 1)
117 126004 : msphib = SIZE(sphi_b, 1)
118 126004 : msphic = SIZE(sphi_c, 1)
119 :
120 126004 : m1 = SIZE(sabc, 1)
121 126004 : m2 = SIZE(sabc, 2)
122 126004 : m3 = SIZE(sabc, 3)
123 126004 : mx = MAX(m2, nsgfb)
124 :
125 : ! ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, nsgfb, m3))
126 630020 : ALLOCATE (tmp(nsgfa, mx, m3 + 1))
127 :
128 126004 : CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, tmp(:, :, 2), nsgfa)
129 1271092 : DO i = 1, m3
130 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, tmp(:, :, i + 1), nsgfa, sphi_b, msphib, &
131 1271092 : 0.0_dp, tmp(:, :, i), nsgfa)
132 : END DO
133 : CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, 1._dp, tmp, nsgfa*mx, sphi_c, msphic, 0.0_dp, &
134 1155536 : abcint, nsgfa*nsgfb)
135 :
136 126004 : DEALLOCATE (tmp)
137 :
138 126004 : CALL timestop(handle)
139 :
140 126004 : END SUBROUTINE abc_contract
141 :
142 : ! **************************************************************************************************
143 : !> \brief contract four-center overlap integrals (a,b,c,d) and transfer
144 : !> to spherical Gaussians
145 : !> \param abcdint contracted, normalized integrals of spherical Gaussians
146 : !> \param sabcd uncontracted, unnormalized integrals of primitive Cartesian Gaussians
147 : !> \param sphi_a contraction matrix for center a
148 : !> \param sphi_b contraction matrix for center b
149 : !> \param sphi_c contraction matrix for center c
150 : !> \param sphi_d contraction matrix for center d
151 : !> \param ncoa number of cartesian orbitals on a
152 : !> \param ncob number of cartesian orbitals on b
153 : !> \param ncoc number of cartesian orbitals on c
154 : !> \param ncod number of cartesian orbitals on d
155 : !> \param nsgfa number of spherical Gaussian functions on a
156 : !> \param nsgfb number of spherical Gaussian functions on b
157 : !> \param nsgfc number of spherical Gaussian functions on c
158 : !> \param nsgfd number of spherical Gaussian functions on d
159 : ! **************************************************************************************************
160 9 : SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, &
161 : ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
162 :
163 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
164 : INTENT(INOUT) :: abcdint
165 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sabcd
166 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sphi_a, sphi_b, sphi_c, sphi_d
167 : INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, &
168 : nsgfc, nsgfd
169 :
170 : CHARACTER(LEN=*), PARAMETER :: routineN = 'abcd_contract'
171 :
172 : INTEGER :: handle, isgfc, isgfd, m1, m2, m3, m4, &
173 : msphia, msphib, msphic, msphid
174 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: temp_cccc, work_cpcc
175 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: temp_cpcc, work_cppc
176 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: cpcc, cppc, cppp
177 :
178 9 : CALL timeset(routineN, handle)
179 :
180 9 : msphia = SIZE(sphi_a, 1)
181 9 : msphib = SIZE(sphi_b, 1)
182 9 : msphic = SIZE(sphi_c, 1)
183 9 : msphid = SIZE(sphi_d, 1)
184 :
185 9 : m1 = SIZE(sabcd, 1)
186 9 : m2 = SIZE(sabcd, 2)
187 9 : m3 = SIZE(sabcd, 3)
188 9 : m4 = SIZE(sabcd, 4)
189 :
190 : ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
191 144 : cpcc(nsgfa, m2, nsgfc, nsgfd))
192 :
193 81 : ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
194 9 : work_cppc = 0._dp
195 9 : temp_cpcc = 0._dp
196 :
197 63 : ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
198 9 : work_cpcc = 0._dp
199 9 : temp_cccc = 0._dp
200 :
201 : CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
202 9 : 0.0_dp, cppp, nsgfa)
203 : CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
204 9 : sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
205 :
206 45 : DO isgfd = 1, nsgfd
207 142164 : work_cppc(:, :, :) = cppc(:, :, :, isgfd)
208 : CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
209 36 : sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
210 20340 : cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
211 189 : DO isgfc = 1, nsgfc
212 20304 : work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
213 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
214 144 : msphib, 0.0_dp, temp_cccc, nsgfa)
215 3060 : abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
216 : END DO
217 : END DO
218 :
219 9 : DEALLOCATE (cpcc, cppc, cppp)
220 9 : DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
221 :
222 9 : CALL timestop(handle)
223 :
224 9 : END SUBROUTINE abcd_contract
225 :
226 : ! **************************************************************************************************
227 : !> \brief 3-center contraction routine from primitive cartesian Gaussians to spherical Gaussian
228 : !> functions using BLAS dgemm.
229 : !> Requires pre-transposition of the sphi_a array. The call-side shall DEALLOCATE buffers
230 : !> end of scope or after last use. This function ALLOCATEs or grows the work buffers
231 : !> as necessary.
232 : !> \param abcint contracted integrals
233 : !> \param sabc uncontracted integrals
234 : !> \param sphi_a assumed to have dimensions nsgfa x ncoa
235 : !> \param sphi_b assumed to have dimensions ncob x nsgfb
236 : !> \param sphi_c assumed to have dimensions ncoc x nsgfc
237 : !> \param ncoa ...
238 : !> \param ncob ...
239 : !> \param ncoc ...
240 : !> \param nsgfa ...
241 : !> \param nsgfb ...
242 : !> \param nsgfc ...
243 : !> \param cpp_buffer Buffer used for intermediate results (automatically allocated).
244 : !> \param ccp_buffer Buffer used for intermediate results (automatically allocated).
245 : !> \param prefac Prefactor which is finally multiplied into abcint (default: 1.0).
246 : !> \param pstfac Factor used to consider initial abcint (default: 0.0).
247 : ! **************************************************************************************************
248 2202416 : SUBROUTINE abc_contract_xsmm(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
249 : nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer, prefac, pstfac)
250 :
251 : REAL(KIND=dp), DIMENSION(:, :, :) :: abcint
252 : REAL(KIND=dp), DIMENSION(*), INTENT(IN), TARGET :: sabc
253 : REAL(KIND=dp), DIMENSION(:, :), CONTIGUOUS, INTENT(IN), TARGET :: sphi_a, sphi_b, sphi_c
254 : INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
255 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, TARGET :: cpp_buffer, ccp_buffer
256 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: prefac, pstfac
257 :
258 : CHARACTER(LEN=*), PARAMETER :: routineN = 'abc_contract_xsmm'
259 :
260 : REAL(KIND=dp) :: alpha, beta
261 : INTEGER(KIND=int_8) :: cpp_size, ccp_size
262 : INTEGER :: handle, i
263 : LOGICAL :: ab_first
264 : #if defined(__LIBXS)
265 : TYPE(libxs_gemm_config_t) :: cfg1, cfg2
266 : INTEGER :: rc1, rc2
267 : #endif
268 :
269 2202416 : CALL timeset(routineN, handle)
270 :
271 2202416 : alpha = 1.0_dp
272 2202416 : IF (PRESENT(prefac)) alpha = prefac
273 :
274 2202416 : beta = 0.0_dp
275 2202416 : IF (PRESENT(pstfac)) beta = pstfac
276 :
277 : ! M*N*K FLOPS are used to decide if contracting (AB)C vs A(BC)
278 2202416 : IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (ncoa*nsgfb*(ncob + nsgfa))) THEN
279 2202416 : cpp_size = nsgfa*ncob
280 2202416 : ab_first = .TRUE.
281 : ELSE
282 655098 : cpp_size = ncoa*nsgfb
283 655098 : ab_first = .FALSE.
284 : END IF
285 :
286 2202416 : ccp_size = nsgfa*nsgfb*ncoc
287 2202416 : IF (.NOT. ALLOCATED(ccp_buffer)) THEN
288 7524 : ALLOCATE (ccp_buffer(ccp_size))
289 2199908 : ELSE IF (SIZE(ccp_buffer) < ccp_size) THEN
290 6542 : DEALLOCATE (ccp_buffer)
291 19626 : ALLOCATE (ccp_buffer(ccp_size))
292 : END IF
293 :
294 2202416 : IF (.NOT. ALLOCATED(cpp_buffer)) THEN
295 7524 : ALLOCATE (cpp_buffer(cpp_size))
296 2199908 : ELSE IF (SIZE(cpp_buffer) < cpp_size) THEN
297 3327 : DEALLOCATE (cpp_buffer)
298 9981 : ALLOCATE (cpp_buffer(cpp_size))
299 : END IF
300 :
301 : #if defined(__LIBXS)
302 2202416 : IF (ab_first) THEN
303 : rc1 = libxs_gemm_dispatch(cfg1, LIBXS_DATATYPE_F64, 'N', 'N', &
304 1547318 : nsgfa, ncob, ncoa, nsgfa, ncoa, nsgfa, 1.0D0, 0.0D0)
305 : rc2 = libxs_gemm_dispatch(cfg2, LIBXS_DATATYPE_F64, 'N', 'N', &
306 1547318 : nsgfa, nsgfb, ncob, nsgfa, ncob, nsgfa, 1.0D0, 0.0D0)
307 : ELSE
308 : rc1 = libxs_gemm_dispatch(cfg1, LIBXS_DATATYPE_F64, 'N', 'N', &
309 655098 : ncoa, nsgfb, ncob, ncoa, ncob, ncoa, 1.0D0, 0.0D0)
310 : rc2 = libxs_gemm_dispatch(cfg2, LIBXS_DATATYPE_F64, 'N', 'N', &
311 655098 : nsgfa, nsgfb, ncoa, nsgfa, ncoa, nsgfa, 1.0D0, 0.0D0)
312 : END IF
313 2202416 : IF (0 /= rc1 .AND. 0 /= rc2) THEN
314 2202416 : IF (ab_first) THEN
315 28105058 : DO i = 0, ncoc - 1
316 : CALL libxs_gemm_call(cfg1, C_LOC(sphi_a(1, 1)), &
317 26557740 : C_LOC(sabc(i*ncoa*ncob + 1)), C_LOC(cpp_buffer(1)))
318 : CALL libxs_gemm_call(cfg2, C_LOC(cpp_buffer(1)), &
319 28105058 : C_LOC(sphi_b(1, 1)), C_LOC(ccp_buffer(i*nsgfa*nsgfb + 1)))
320 : END DO
321 : ELSE
322 11375938 : DO i = 0, ncoc - 1
323 : CALL libxs_gemm_call(cfg1, C_LOC(sabc(i*ncoa*ncob + 1)), &
324 10720840 : C_LOC(sphi_b(1, 1)), C_LOC(cpp_buffer(1)))
325 : CALL libxs_gemm_call(cfg2, C_LOC(sphi_a(1, 1)), &
326 11375938 : C_LOC(cpp_buffer(1)), C_LOC(ccp_buffer(i*nsgfa*nsgfb + 1)))
327 : END DO
328 : END IF
329 : ELSE
330 : #endif
331 0 : IF (ab_first) THEN ! (AB)C
332 0 : DO i = 0, ncoc - 1
333 : CALL dgemm("N", "N", nsgfa, ncob, ncoa, 1.0_dp, sphi_a, nsgfa, sabc(i*ncoa*ncob + 1), &
334 0 : ncoa, 0.0_dp, cpp_buffer, nsgfa)
335 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1.0_dp, cpp_buffer, nsgfa, sphi_b, ncob, 0.0_dp, &
336 0 : ccp_buffer(i*nsgfa*nsgfb + 1), nsgfa)
337 : END DO
338 : ELSE ! A(BC)
339 0 : DO i = 0, ncoc - 1
340 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sabc(i*ncoa*ncob + 1), ncoa, sphi_b, &
341 0 : ncob, 0.0_dp, cpp_buffer, ncoa)
342 : CALL dgemm("N", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a, nsgfa, cpp_buffer, ncoa, 0.0_dp, &
343 0 : ccp_buffer(i*nsgfa*nsgfb + 1), nsgfa)
344 : END DO
345 : END IF
346 : #if defined(__LIBXS)
347 : END IF
348 : #endif
349 : ! contractions over c: [nsgfa*nsgfb,ncoc] x [ncoc,nsgfc] -> [sgfa*nsgfb,nsgfc]
350 : CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, alpha, ccp_buffer, nsgfa*nsgfb, &
351 2202416 : sphi_c, ncoc, beta, abcint, nsgfa*nsgfb)
352 :
353 2202416 : CALL timestop(handle)
354 :
355 2202416 : END SUBROUTINE abc_contract_xsmm
356 :
357 : END MODULE ai_contraction_sphi
|
