(git:374b731)
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-2024 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
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. Exploits LIBXSMM for performance, falls back to BLAS if LIBXSMM is not available.
234!> Requires pre-allocation of work buffers and pre-transposition of the sphi_a array. Requires
235!> the LIBXSMM library to be initialized somewhere before this routine is called.
236!> \param abcint contracted integrals
237!> \param sabc uncontracted integrals
238!> \param sphi_a assumed to have dimensions nsgfa x ncoa
239!> \param sphi_b assumed to have dimensions ncob x nsgfb
240!> \param sphi_c assumed to have dimensions ncoc x nsgfc
241!> \param ncoa ...
242!> \param ncob ...
243!> \param ncoc ...
244!> \param nsgfa ...
245!> \param nsgfb ...
246!> \param nsgfc ...
247!> \param cpp_buffer Temporary buffer used for intermediate results.
248!> \param ccp_buffer Temporary buffer used for intermediate results.
249!> \param prefac Prefactor which is finally multiplied into abcint.
250!> \note tested from version 1.9.0 of libxsmm
251! **************************************************************************************************
252 SUBROUTINE abc_contract_xsmm(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
253 nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer, prefac)
254
255 REAL(kind=dp), DIMENSION(*) :: abcint
256 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: sabc
257 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sphi_a, sphi_b, sphi_c
258 INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
259 REAL(kind=dp), DIMENSION(nsgfa*ncob), INTENT(OUT) :: cpp_buffer
260 REAL(kind=dp), DIMENSION(nsgfa*nsgfb*ncoc), INTENT(OUT) :: ccp_buffer
261 REAL(kind=dp), INTENT(IN), OPTIONAL :: prefac
262
263 CHARACTER(LEN=*), PARAMETER :: routinen = 'abc_contract_xsmm'
264
265 REAL(kind=dp) :: alpha
266 INTEGER :: handle, i
267 LOGICAL :: ab
268#if defined(__LIBXSMM)
269 TYPE(libxsmm_dmmfunction) :: xmm1, xmm2
270#endif
271
272 ! M*N*K FLOPS can be used to decide if contracting (AB)C vs A(BC)
273 ! ab = (nsgfa*ncob*(ncoa + nsgfb)) <= (ncoa*nsgfb*(ncob + nsgfa))
274 ! however, current interface dictates size of cpp_buffer
275 ab = (nsgfa*ncob) < (ncoa*nsgfb) .OR. (ncoa + nsgfb) <= (ncob + nsgfa)
276
277 CALL timeset(routinen, handle)
278
279 alpha = 1.0_dp
280 IF (PRESENT(prefac)) alpha = prefac
281
282 ! loop over the last index of the matrix and call LIBXSMM/BLAS to contract over a and b
283#if defined(__LIBXSMM)
284 IF (ab) THEN ! (AB)C: dispatch kernels
285 CALL libxsmm_dispatch(xmm1, nsgfa, ncob, ncoa, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
286 CALL libxsmm_dispatch(xmm2, nsgfa, nsgfb, ncob, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
287 ELSE ! A(BC): dispatch kernels
288 CALL libxsmm_dispatch(xmm1, ncoa, nsgfb, ncob, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
289 CALL libxsmm_dispatch(xmm2, nsgfa, nsgfb, ncoa, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
290 END IF
291
292 IF (libxsmm_available(xmm1) .AND. libxsmm_available(xmm2)) THEN
293 IF (ab) THEN ! (AB)C
294 DO i = 1, ncoc ! contractions over a and b
295 ! [nsgfa,ncoa] x [ncoa,ncob] -> [nsgfa,ncob]
296 CALL libxsmm_dmmcall(xmm1, sphi_a, sabc((i - 1)*ncoa*ncob + 1), cpp_buffer)
297 ! [nsgfa,ncob] x [ncob,nsgfb] -> [nsgfa,nsgfb]
298 CALL libxsmm_dmmcall(xmm2, cpp_buffer, sphi_b, ccp_buffer((i - 1)*nsgfa*nsgfb + 1))
299 END DO
300 ELSE ! A(BC)
301 DO i = 1, ncoc ! contractions over a and b
302 ! [ncoa,ncob] x [ncob,nsgfb] -> [ncoa,nsgfb]
303 CALL libxsmm_dmmcall(xmm1, sabc((i - 1)*ncoa*ncob + 1), sphi_b, cpp_buffer)
304 ! [nsgfa,ncoa] x [ncoa,nsgfb] -> [nsgfa,nsgfb]
305 CALL libxsmm_dmmcall(xmm2, sphi_a, cpp_buffer, ccp_buffer((i - 1)*nsgfa*nsgfb + 1))
306 END DO
307 END IF
308 ELSE
309#endif
310 IF (ab) THEN ! (AB)C
311 DO i = 1, ncoc ! contractions over a and b
312 CALL dgemm("N", "N", nsgfa, ncob, ncoa, 1.0_dp, sphi_a, nsgfa, sabc((i - 1)*ncoa*ncob + 1), &
313 ncoa, 0.0_dp, cpp_buffer, nsgfa)
314 CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1.0_dp, cpp_buffer, nsgfa, sphi_b, ncob, 0.0_dp, &
315 ccp_buffer((i - 1)*nsgfa*nsgfb + 1), nsgfa)
316 END DO
317 ELSE ! A(BC)
318 DO i = 1, ncoc ! contractions over a and b
319 CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sabc((i - 1)*ncoa*ncob + 1), ncoa, sphi_b, &
320 ncob, 0.0_dp, cpp_buffer, ncoa)
321 CALL dgemm("N", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a, nsgfa, cpp_buffer, ncoa, 0.0_dp, &
322 ccp_buffer((i - 1)*nsgfa*nsgfb + 1), nsgfa)
323 END DO
324 END IF
325#if defined(__LIBXSMM)
326 END IF
327#endif
328 ! last contraction, over c (larger MM using BLAS)
329 CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, alpha, ccp_buffer, nsgfa*nsgfb, sphi_c, ncoc, &
330 0.0_dp, abcint, nsgfa*nsgfb)
331
332 CALL timestop(handle)
333
334 END SUBROUTINE abc_contract_xsmm
335
336END 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:195
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_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)
3-center contraction routine from primitive cartesian Gaussians to spherical Gaussian functions....
Definition ai_contraction_sphi.F:254
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::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::dp
integer, parameter, public dp
Definition kinds.F:34