Line data Source code
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 : ! **************************************************************************************************
15 : MODULE 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 :
38 : CONTAINS
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 1108678 : 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 1108678 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cpp
61 :
62 1108678 : CALL timeset(routineN, handle)
63 :
64 1108678 : CPASSERT(ncob <= SIZE(sab, 2))
65 :
66 1108678 : IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (nsgfb*ncoa*(ncob + nsgfa))) THEN ! (sphi_a x sab) x sphi_b
67 3028304 : ALLOCATE (cpp(nsgfa, ncob))
68 : ! [nsgfa,ncoa] x [ncoa,ncob] -> [nsgfa,ncob]
69 22289401 : 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 52976896 : abint, SIZE(abint, 1))
73 : ELSE ! sphi_a x (sab x sphi_b)
74 1406408 : ALLOCATE (cpp(ncoa, nsgfb))
75 : ! [ncoa,ncob] x [ncob,nsgfb] -> [ncoa,nsgfb]
76 7604846 : 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 13521000 : abint, SIZE(abint, 1))
80 : END IF
81 :
82 1108678 : DEALLOCATE (cpp)
83 :
84 1108678 : CALL timestop(handle)
85 :
86 1108678 : 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 96316 : 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 96316 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: tmp
115 :
116 96316 : CALL timeset(routineN, handle)
117 :
118 96316 : CPASSERT(SIZE(abcint, 1) == nsgfa)
119 96316 : CPASSERT(SIZE(abcint, 2) == nsgfb)
120 :
121 96316 : msphia = SIZE(sphi_a, 1)
122 96316 : msphib = SIZE(sphi_b, 1)
123 96316 : msphic = SIZE(sphi_c, 1)
124 :
125 96316 : m1 = SIZE(sabc, 1)
126 96316 : m2 = SIZE(sabc, 2)
127 96316 : m3 = SIZE(sabc, 3)
128 96316 : mx = MAX(m2, nsgfb)
129 :
130 : ! ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, nsgfb, m3))
131 481580 : ALLOCATE (tmp(nsgfa, mx, m3 + 1))
132 :
133 96316 : CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, tmp(:, :, 2), nsgfa)
134 1076476 : DO i = 1, m3
135 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, tmp(:, :, i + 1), nsgfa, sphi_b, msphib, &
136 1076476 : 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 1125848 : abcint, nsgfa*nsgfb)
140 :
141 96316 : DEALLOCATE (tmp)
142 :
143 96316 : CALL timestop(handle)
144 :
145 96316 : 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 9 : 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 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: temp_cccc, work_cpcc
180 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: temp_cpcc, work_cppc
181 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: cpcc, cppc, cppp
182 :
183 9 : CALL timeset(routineN, handle)
184 :
185 9 : msphia = SIZE(sphi_a, 1)
186 9 : msphib = SIZE(sphi_b, 1)
187 9 : msphic = SIZE(sphi_c, 1)
188 9 : msphid = SIZE(sphi_d, 1)
189 :
190 9 : m1 = SIZE(sabcd, 1)
191 9 : m2 = SIZE(sabcd, 2)
192 9 : m3 = SIZE(sabcd, 3)
193 9 : m4 = SIZE(sabcd, 4)
194 :
195 : ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
196 144 : cpcc(nsgfa, m2, nsgfc, nsgfd))
197 :
198 81 : ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
199 35541 : work_cppc = 0._dp
200 5085 : temp_cpcc = 0._dp
201 :
202 63 : ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
203 1269 : work_cpcc = 0._dp
204 189 : temp_cccc = 0._dp
205 :
206 : CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
207 9 : 0.0_dp, cppp, nsgfa)
208 : CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
209 9 : sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
210 :
211 45 : DO isgfd = 1, nsgfd
212 142164 : work_cppc(:, :, :) = cppc(:, :, :, isgfd)
213 : CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
214 36 : sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
215 20340 : cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
216 189 : DO isgfc = 1, nsgfc
217 20304 : work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
218 : CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
219 144 : msphib, 0.0_dp, temp_cccc, nsgfa)
220 3060 : abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
221 : END DO
222 : END DO
223 :
224 9 : DEALLOCATE (cpcc, cppc, cppp)
225 9 : DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
226 :
227 9 : CALL timestop(handle)
228 :
229 9 : 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 2480509 : 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 2480509 : CALL timeset(routineN, handle)
274 :
275 2480509 : alpha = 1.0_dp
276 2480509 : IF (PRESENT(prefac)) alpha = prefac
277 :
278 2480509 : beta = 0.0_dp
279 2480509 : IF (PRESENT(pstfac)) beta = pstfac
280 :
281 : ! M*N*K FLOPS are used to decide if contracting (AB)C vs A(BC)
282 2480509 : IF ((nsgfa*ncob*(ncoa + nsgfb)) <= (ncoa*nsgfb*(ncob + nsgfa))) THEN
283 1727787 : cpp_size = nsgfa*ncob
284 1727787 : ab_first = .TRUE.
285 : ELSE
286 752722 : cpp_size = ncoa*nsgfb
287 752722 : ab_first = .FALSE.
288 : END IF
289 :
290 2480509 : ccp_size = nsgfa*nsgfb*ncoc
291 2480509 : IF (.NOT. ALLOCATED(ccp_buffer)) THEN
292 5322 : ALLOCATE (ccp_buffer(ccp_size))
293 2478735 : ELSE IF (SIZE(ccp_buffer) < ccp_size) THEN
294 6884 : DEALLOCATE (ccp_buffer)
295 20652 : ALLOCATE (ccp_buffer(ccp_size))
296 : END IF
297 :
298 2480509 : IF (.NOT. ALLOCATED(cpp_buffer)) THEN
299 5322 : ALLOCATE (cpp_buffer(cpp_size))
300 2478735 : ELSE IF (SIZE(cpp_buffer) < cpp_size) THEN
301 4525 : DEALLOCATE (cpp_buffer)
302 13575 : 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 2480509 : IF (ab_first) THEN ! (AB)C: dispatch kernels
308 1727787 : CALL libxsmm_dispatch(xmm1, nsgfa, ncob, ncoa, beta=0.0_dp, prefetch=LIBXSMM_PREFETCH_NONE)
309 1727787 : CALL libxsmm_dispatch(xmm2, nsgfa, nsgfb, ncob, beta=0.0_dp, prefetch=LIBXSMM_PREFETCH_NONE)
310 : ELSE ! A(BC): dispatch kernels
311 752722 : CALL libxsmm_dispatch(xmm1, ncoa, nsgfb, ncob, beta=0.0_dp, prefetch=LIBXSMM_PREFETCH_NONE)
312 752722 : CALL libxsmm_dispatch(xmm2, nsgfa, nsgfb, ncoa, beta=0.0_dp, prefetch=LIBXSMM_PREFETCH_NONE)
313 : END IF
314 :
315 2480509 : IF (libxsmm_available(xmm1) .AND. libxsmm_available(xmm2)) THEN
316 2480509 : IF (ab_first) THEN ! (AB)C
317 25911975 : DO i = 0, ncoc - 1 ! contractions over a and b
318 : ! [nsgfa,ncoa] x [ncoa,ncob] -> [nsgfa,ncob]
319 24184188 : CALL libxsmm_dmmcall(xmm1, sphi_a, sabc(i*ncoa*ncob + 1), cpp_buffer)
320 : ! [nsgfa,ncob] x [ncob,nsgfb] -> [nsgfa,nsgfb]
321 25911975 : CALL libxsmm_dmmcall(xmm2, cpp_buffer, sphi_b, ccp_buffer(i*nsgfa*nsgfb + 1))
322 : END DO
323 : ELSE ! A(BC)
324 10425276 : DO i = 0, ncoc - 1 ! contractions over a and b
325 : ! [ncoa,ncob] x [ncob,nsgfb] -> [ncoa,nsgfb]
326 9672554 : CALL libxsmm_dmmcall(xmm1, sabc(i*ncoa*ncob + 1), sphi_b, cpp_buffer)
327 : ! [nsgfa,ncoa] x [ncoa,nsgfb] -> [nsgfa,nsgfb]
328 10425276 : 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 0 : IF (ab_first) THEN ! (AB)C
334 0 : 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 0 : 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 0 : ccp_buffer(i*nsgfa*nsgfb + 1), nsgfa) ! [nsgfa,ncob] x [ncob,nsgfb] -> [nsgfa,nsgfb]
339 : END DO
340 : ELSE ! A(BC)
341 0 : 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 0 : 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 0 : 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 2480509 : sphi_c, ncoc, beta, abcint, nsgfa*nsgfb)
354 :
355 2480509 : CALL timestop(handle)
356 :
357 2480509 : END SUBROUTINE abc_contract_xsmm
358 :
359 : END MODULE ai_contraction_sphi
|
