(git:3c69fd1)
Loading...
Searching...
No Matches
ai_contraction_sphi.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief 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! **************************************************************************************************
15MODULE 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
33CONTAINS
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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cpp
56
57 CALL timeset(routinen, handle)
58
59 cpassert(ncob <= SIZE(sab, 2))
60
61 IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (nsgfb*ncoa*(ncob + nsgfa))) THEN ! (sphi_a x sab) x sphi_b
62 ALLOCATE (cpp(nsgfa, ncob))
63 ! [nsgfa,ncoa] x [ncoa,ncob] -> [nsgfa,ncob]
64 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 abint, SIZE(abint, 1))
68 ELSE ! sphi_a x (sab x sphi_b)
69 ALLOCATE (cpp(ncoa, nsgfb))
70 ! [ncoa,ncob] x [ncob,nsgfb] -> [ncoa,nsgfb]
71 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 abint, SIZE(abint, 1))
75 END IF
76
77 DEALLOCATE (cpp)
78
79 CALL timestop(handle)
80
81 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: tmp
110
111 CALL timeset(routinen, handle)
112
113 cpassert(SIZE(abcint, 1) == nsgfa)
114 cpassert(SIZE(abcint, 2) == nsgfb)
115
116 msphia = SIZE(sphi_a, 1)
117 msphib = SIZE(sphi_b, 1)
118 msphic = SIZE(sphi_c, 1)
119
120 m1 = SIZE(sabc, 1)
121 m2 = SIZE(sabc, 2)
122 m3 = SIZE(sabc, 3)
123 mx = max(m2, nsgfb)
124
125 ! ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, nsgfb, m3))
126 ALLOCATE (tmp(nsgfa, mx, m3 + 1))
127
128 CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, tmp(:, :, 2), nsgfa)
129 DO i = 1, m3
130 CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, tmp(:, :, i + 1), nsgfa, sphi_b, msphib, &
131 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 abcint, nsgfa*nsgfb)
135
136 DEALLOCATE (tmp)
137
138 CALL timestop(handle)
139
140 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: temp_cccc, work_cpcc
175 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: temp_cpcc, work_cppc
176 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: cpcc, cppc, cppp
177
178 CALL timeset(routinen, handle)
179
180 msphia = SIZE(sphi_a, 1)
181 msphib = SIZE(sphi_b, 1)
182 msphic = SIZE(sphi_c, 1)
183 msphid = SIZE(sphi_d, 1)
184
185 m1 = SIZE(sabcd, 1)
186 m2 = SIZE(sabcd, 2)
187 m3 = SIZE(sabcd, 3)
188 m4 = SIZE(sabcd, 4)
189
190 ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
191 cpcc(nsgfa, m2, nsgfc, nsgfd))
192
193 ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
194 work_cppc = 0._dp
195 temp_cpcc = 0._dp
196
197 ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
198 work_cpcc = 0._dp
199 temp_cccc = 0._dp
200
201 CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
202 0.0_dp, cppp, nsgfa)
203 CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
204 sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
205
206 DO isgfd = 1, nsgfd
207 work_cppc(:, :, :) = cppc(:, :, :, isgfd)
208 CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
209 sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
210 cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
211 DO isgfc = 1, nsgfc
212 work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
213 CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
214 msphib, 0.0_dp, temp_cccc, nsgfa)
215 abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
216 END DO
217 END DO
218
219 DEALLOCATE (cpcc, cppc, cppp)
220 DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
221
222 CALL timestop(handle)
223
224 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 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 CALL timeset(routinen, handle)
270
271 alpha = 1.0_dp
272 IF (PRESENT(prefac)) alpha = prefac
273
274 beta = 0.0_dp
275 IF (PRESENT(pstfac)) beta = pstfac
276
277 ! M*N*K FLOPS are used to decide if contracting (AB)C vs A(BC)
278 IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (ncoa*nsgfb*(ncob + nsgfa))) THEN
279 cpp_size = nsgfa*ncob
280 ab_first = .true.
281 ELSE
282 cpp_size = ncoa*nsgfb
283 ab_first = .false.
284 END IF
285
286 ccp_size = nsgfa*nsgfb*ncoc
287 IF (.NOT. ALLOCATED(ccp_buffer)) THEN
288 ALLOCATE (ccp_buffer(ccp_size))
289 ELSE IF (SIZE(ccp_buffer) < ccp_size) THEN
290 DEALLOCATE (ccp_buffer)
291 ALLOCATE (ccp_buffer(ccp_size))
292 END IF
293
294 IF (.NOT. ALLOCATED(cpp_buffer)) THEN
295 ALLOCATE (cpp_buffer(cpp_size))
296 ELSE IF (SIZE(cpp_buffer) < cpp_size) THEN
297 DEALLOCATE (cpp_buffer)
298 ALLOCATE (cpp_buffer(cpp_size))
299 END IF
300
301#if defined(__LIBXS)
302 IF (ab_first) THEN
303 rc1 = libxs_gemm_dispatch(cfg1, libxs_datatype_f64, 'N', 'N', &
304 nsgfa, ncob, ncoa, nsgfa, ncoa, nsgfa, 1.0d0, 0.0d0)
305 rc2 = libxs_gemm_dispatch(cfg2, libxs_datatype_f64, 'N', 'N', &
306 nsgfa, nsgfb, ncob, nsgfa, ncob, nsgfa, 1.0d0, 0.0d0)
307 ELSE
308 rc1 = libxs_gemm_dispatch(cfg1, libxs_datatype_f64, 'N', 'N', &
309 ncoa, nsgfb, ncob, ncoa, ncob, ncoa, 1.0d0, 0.0d0)
310 rc2 = libxs_gemm_dispatch(cfg2, libxs_datatype_f64, 'N', 'N', &
311 nsgfa, nsgfb, ncoa, nsgfa, ncoa, nsgfa, 1.0d0, 0.0d0)
312 END IF
313 IF (0 /= rc1 .AND. 0 /= rc2) THEN
314 IF (ab_first) THEN
315 DO i = 0, ncoc - 1
316 CALL libxs_gemm_call(cfg1, c_loc(sphi_a(1, 1)), &
317 c_loc(sabc(i*ncoa*ncob + 1)), c_loc(cpp_buffer(1)))
318 CALL libxs_gemm_call(cfg2, c_loc(cpp_buffer(1)), &
319 c_loc(sphi_b(1, 1)), c_loc(ccp_buffer(i*nsgfa*nsgfb + 1)))
320 END DO
321 ELSE
322 DO i = 0, ncoc - 1
323 CALL libxs_gemm_call(cfg1, c_loc(sabc(i*ncoa*ncob + 1)), &
324 c_loc(sphi_b(1, 1)), c_loc(cpp_buffer(1)))
325 CALL libxs_gemm_call(cfg2, c_loc(sphi_a(1, 1)), &
326 c_loc(cpp_buffer(1)), c_loc(ccp_buffer(i*nsgfa*nsgfb + 1)))
327 END DO
328 END IF
329 ELSE
330#endif
331 IF (ab_first) THEN ! (AB)C
332 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 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 ccp_buffer(i*nsgfa*nsgfb + 1), nsgfa)
337 END DO
338 ELSE ! A(BC)
339 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 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 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 sphi_c, ncoc, beta, abcint, nsgfa*nsgfb)
352
353 CALL timestop(handle)
354
355 END SUBROUTINE abc_contract_xsmm
356
357END MODULE ai_contraction_sphi
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
ai_contraction_sphi
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
Definition ai_contraction_sphi.F:15
ai_contraction_sphi::abc_contract
subroutine, public abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
contract three-center overlap integrals (a,b,c) and transfer to spherical Gaussians
Definition ai_contraction_sphi.F:100
ai_contraction_sphi::abcd_contract
subroutine, public abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
contract four-center overlap integrals (a,b,c,d) and transfer to spherical Gaussians
Definition ai_contraction_sphi.F:162
ai_contraction_sphi::abc_contract_xsmm
subroutine, public abc_contract_xsmm(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer, prefac, pstfac)
3-center contraction routine from primitive cartesian Gaussians to spherical Gaussian functions using...
Definition ai_contraction_sphi.F:250
ai_contraction_sphi::ab_contract
subroutine, public ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
contract overlap integrals (a,b) and transfer to spherical Gaussians
Definition ai_contraction_sphi.F:47
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34