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