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 1292952 : 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 1292952 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cpp
56 :
57 1292952 : CALL timeset(routineN, handle)
58 :
59 1292952 : CPASSERT(ncob <= SIZE(sab, 2))
60 :
61 1292952 : IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (nsgfb*ncoa*(ncob + nsgfa))) THEN ! (sphi_a x sab) x sphi_b
62 3515364 : ALLOCATE (cpp(nsgfa, ncob))
63 : ! [nsgfa,ncoa] x [ncoa,ncob] -> [nsgfa,ncob]
64 24712188 : 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 55412391 : abint, SIZE(abint, 1))
68 : ELSE ! sphi_a x (sab x sphi_b)
69 1656444 : ALLOCATE (cpp(ncoa, nsgfb))
70 : ! [ncoa,ncob] x [ncob,nsgfb] -> [ncoa,nsgfb]
71 8345071 : 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 14508899 : abint, SIZE(abint, 1))
75 : END IF
76 :
77 1292952 : DEALLOCATE (cpp)
78 :
79 1292952 : CALL timestop(handle)
80 :
81 1292952 : 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 207072 : 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 207072 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: tmp
110 :
111 207072 : CALL timeset(routineN, handle)
112 :
113 207072 : CPASSERT(SIZE(abcint, 1) == nsgfa)
114 207072 : CPASSERT(SIZE(abcint, 2) == nsgfb)
115 :
116 207072 : msphia = SIZE(sphi_a, 1)
117 207072 : msphib = SIZE(sphi_b, 1)
118 207072 : msphic = SIZE(sphi_c, 1)
119 :
120 207072 : m1 = SIZE(sabc, 1)
121 207072 : m2 = SIZE(sabc, 2)
122 207072 : m3 = SIZE(sabc, 3)
123 207072 : mx = MAX(m2, nsgfb)
124 :
125 : ! ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, nsgfb, m3))
126 1035360 : ALLOCATE (tmp(nsgfa, mx, m3 + 1))
127 :
128 207072 : CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, tmp(:, :, 2), nsgfa)
129 1959016 : DO i = 1, m3
130 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, tmp(:, :, i + 1), nsgfa, sphi_b, msphib, &
131 1959016 : 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 1236604 : abcint, nsgfa*nsgfb)
135 :
136 207072 : DEALLOCATE (tmp)
137 :
138 207072 : CALL timestop(handle)
139 :
140 207072 : 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
175 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpcc
176 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: 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 162 : cpcc(nsgfa, m2, nsgfc), temp_cccc(nsgfa, nsgfb))
192 :
193 : CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
194 9 : 0.0_dp, cppp, nsgfa)
195 : CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
196 9 : sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
197 :
198 45 : DO isgfd = 1, nsgfd
199 : ! Full leading-dimension slices are contiguous; reuse cpcc for each spherical d.
200 : CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, cppc(:, :, :, isgfd), nsgfa*m2, &
201 36 : sphi_c, msphic, 0.0_dp, cpcc, nsgfa*m2)
202 189 : DO isgfc = 1, nsgfc
203 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpcc(:, :, isgfc), nsgfa, sphi_b, &
204 144 : msphib, 0.0_dp, temp_cccc, nsgfa)
205 : ! The caller's output block can be strided, so scatter from a contiguous buffer.
206 3060 : abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
207 : END DO
208 : END DO
209 :
210 9 : DEALLOCATE (cpcc, cppc, cppp, temp_cccc)
211 :
212 9 : CALL timestop(handle)
213 :
214 9 : END SUBROUTINE abcd_contract
215 :
216 : ! **************************************************************************************************
217 : !> \brief 3-center contraction routine from primitive cartesian Gaussians to spherical Gaussian
218 : !> functions using BLAS dgemm.
219 : !> Requires pre-transposition of the sphi_a array. The call-side shall DEALLOCATE buffers
220 : !> end of scope or after last use. This function ALLOCATEs or grows the work buffers
221 : !> as necessary.
222 : !> \param abcint contracted integrals
223 : !> \param sabc uncontracted integrals
224 : !> \param sphi_a assumed to have dimensions nsgfa x ncoa
225 : !> \param sphi_b assumed to have dimensions ncob x nsgfb
226 : !> \param sphi_c assumed to have dimensions ncoc x nsgfc
227 : !> \param ncoa ...
228 : !> \param ncob ...
229 : !> \param ncoc ...
230 : !> \param nsgfa ...
231 : !> \param nsgfb ...
232 : !> \param nsgfc ...
233 : !> \param cpp_buffer Buffer used for intermediate results (automatically allocated).
234 : !> \param ccp_buffer Buffer used for intermediate results (automatically allocated).
235 : !> \param prefac Prefactor which is finally multiplied into abcint (default: 1.0).
236 : !> \param pstfac Factor used to consider initial abcint (default: 0.0).
237 : ! **************************************************************************************************
238 6850397 : SUBROUTINE abc_contract_xsmm(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
239 : nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer, prefac, pstfac)
240 :
241 : REAL(KIND=dp), DIMENSION(:, :, :) :: abcint
242 : REAL(KIND=dp), DIMENSION(*), INTENT(IN), TARGET :: sabc
243 : REAL(KIND=dp), DIMENSION(:, :), CONTIGUOUS, INTENT(IN), TARGET :: sphi_a, sphi_b, sphi_c
244 : INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
245 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, TARGET :: cpp_buffer, ccp_buffer
246 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: prefac, pstfac
247 :
248 : CHARACTER(LEN=*), PARAMETER :: routineN = 'abc_contract_xsmm'
249 :
250 : REAL(KIND=dp) :: alpha, beta
251 : INTEGER(KIND=int_8) :: cpp_size, ccp_size
252 : INTEGER :: handle, i
253 : LOGICAL :: ab_first
254 : #if defined(__LIBXS)
255 : TYPE(libxs_gemm_config_t) :: cfg1, cfg2
256 : INTEGER :: rc1, rc2
257 : #endif
258 :
259 6850397 : CALL timeset(routineN, handle)
260 :
261 6850397 : alpha = 1.0_dp
262 6850397 : IF (PRESENT(prefac)) alpha = prefac
263 :
264 6850397 : beta = 0.0_dp
265 6850397 : IF (PRESENT(pstfac)) beta = pstfac
266 :
267 : ! M*N*K FLOPS are used to decide if contracting (AB)C vs A(BC)
268 6850397 : IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (ncoa*nsgfb*(ncob + nsgfa))) THEN
269 6850397 : cpp_size = nsgfa*ncob
270 6850397 : ab_first = .TRUE.
271 : ELSE
272 1713235 : cpp_size = ncoa*nsgfb
273 1713235 : ab_first = .FALSE.
274 : END IF
275 :
276 6850397 : ccp_size = nsgfa*nsgfb*ncoc
277 6850397 : IF (.NOT. ALLOCATED(ccp_buffer)) THEN
278 26622 : ALLOCATE (ccp_buffer(ccp_size))
279 6841523 : ELSE IF (SIZE(ccp_buffer) < ccp_size) THEN
280 35320 : DEALLOCATE (ccp_buffer)
281 105960 : ALLOCATE (ccp_buffer(ccp_size))
282 : END IF
283 :
284 6850397 : IF (.NOT. ALLOCATED(cpp_buffer)) THEN
285 26622 : ALLOCATE (cpp_buffer(cpp_size))
286 6841523 : ELSE IF (SIZE(cpp_buffer) < cpp_size) THEN
287 12541 : DEALLOCATE (cpp_buffer)
288 37623 : ALLOCATE (cpp_buffer(cpp_size))
289 : END IF
290 :
291 : #if defined(__LIBXS)
292 6850397 : IF (ab_first) THEN
293 : rc1 = libxs_gemm_dispatch(cfg1, LIBXS_DATATYPE_F64, 'N', 'N', &
294 5137162 : nsgfa, ncob, ncoa, nsgfa, ncoa, nsgfa, 1.0D0, 0.0D0)
295 : rc2 = libxs_gemm_dispatch(cfg2, LIBXS_DATATYPE_F64, 'N', 'N', &
296 5137162 : nsgfa, nsgfb, ncob, nsgfa, ncob, nsgfa, 1.0D0, 0.0D0)
297 : ELSE
298 : rc1 = libxs_gemm_dispatch(cfg1, LIBXS_DATATYPE_F64, 'N', 'N', &
299 1713235 : ncoa, nsgfb, ncob, ncoa, ncob, ncoa, 1.0D0, 0.0D0)
300 : rc2 = libxs_gemm_dispatch(cfg2, LIBXS_DATATYPE_F64, 'N', 'N', &
301 1713235 : nsgfa, nsgfb, ncoa, nsgfa, ncoa, nsgfa, 1.0D0, 0.0D0)
302 : END IF
303 6850397 : IF (0 /= rc1 .AND. 0 /= rc2) THEN
304 6850397 : IF (ab_first) THEN
305 45595316 : DO i = 0, ncoc - 1
306 : CALL libxs_gemm_call(cfg1, C_LOC(sphi_a(1, 1)), &
307 40458154 : C_LOC(sabc(i*ncoa*ncob + 1)), C_LOC(cpp_buffer(1)))
308 : CALL libxs_gemm_call(cfg2, C_LOC(cpp_buffer(1)), &
309 45595316 : C_LOC(sphi_b(1, 1)), C_LOC(ccp_buffer(i*nsgfa*nsgfb + 1)))
310 : END DO
311 : ELSE
312 16781137 : DO i = 0, ncoc - 1
313 : CALL libxs_gemm_call(cfg1, C_LOC(sabc(i*ncoa*ncob + 1)), &
314 15067902 : C_LOC(sphi_b(1, 1)), C_LOC(cpp_buffer(1)))
315 : CALL libxs_gemm_call(cfg2, C_LOC(sphi_a(1, 1)), &
316 16781137 : C_LOC(cpp_buffer(1)), C_LOC(ccp_buffer(i*nsgfa*nsgfb + 1)))
317 : END DO
318 : END IF
319 : ELSE
320 : #endif
321 0 : IF (ab_first) THEN ! (AB)C
322 0 : DO i = 0, ncoc - 1
323 : CALL dgemm("N", "N", nsgfa, ncob, ncoa, 1.0_dp, sphi_a, nsgfa, sabc(i*ncoa*ncob + 1), &
324 0 : ncoa, 0.0_dp, cpp_buffer, nsgfa)
325 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1.0_dp, cpp_buffer, nsgfa, sphi_b, ncob, 0.0_dp, &
326 0 : ccp_buffer(i*nsgfa*nsgfb + 1), nsgfa)
327 : END DO
328 : ELSE ! A(BC)
329 0 : DO i = 0, ncoc - 1
330 : CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sabc(i*ncoa*ncob + 1), ncoa, sphi_b, &
331 0 : ncob, 0.0_dp, cpp_buffer, ncoa)
332 : CALL dgemm("N", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a, nsgfa, cpp_buffer, ncoa, 0.0_dp, &
333 0 : ccp_buffer(i*nsgfa*nsgfb + 1), nsgfa)
334 : END DO
335 : END IF
336 : #if defined(__LIBXS)
337 : END IF
338 : #endif
339 : ! contractions over c: [nsgfa*nsgfb,ncoc] x [ncoc,nsgfc] -> [sgfa*nsgfb,nsgfc]
340 : CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, alpha, ccp_buffer, nsgfa*nsgfb, &
341 6850397 : sphi_c, ncoc, beta, abcint, nsgfa*nsgfb)
342 :
343 6850397 : CALL timestop(handle)
344 :
345 6850397 : END SUBROUTINE abc_contract_xsmm
346 :
347 : END MODULE ai_contraction_sphi
|
