(git:0de0cc2)
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 libxsmm_abc_contract routine, A. Bussy (04.2020)
13 !> \author Dorothea Golze (05.2016)
14 ! **************************************************************************************************
16 
17  USE kinds, ONLY: dp
18 #if(__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 
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  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, m1, m2, msphia, msphib, nn
60  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cpp
61 
62  CALL timeset(routinen, handle)
63 
64  msphia = SIZE(sphi_a, 1)
65  msphib = SIZE(sphi_b, 1)
66 
67  m1 = SIZE(sab, 1)
68  m2 = SIZE(sab, 2)
69 
70  nn = SIZE(abint, 1)
71 
72  ALLOCATE (cpp(nsgfa, m2))
73 
74  CALL dgemm("T", "N", nsgfa, m2, ncoa, 1._dp, sphi_a, msphia, sab, m1, 0.0_dp, cpp, nsgfa)
75  CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpp, nsgfa, sphi_b, msphib, 0.0_dp, &
76  abint, nn)
77 
78  DEALLOCATE (cpp)
79 
80  CALL timestop(handle)
81 
82  END SUBROUTINE ab_contract
83 
84 ! **************************************************************************************************
85 !> \brief contract three-center overlap integrals (a,b,c) and transfer
86 !> to spherical Gaussians
87 !> \param abcint contracted, normalized integrals of spherical Gaussians
88 !> \param sabc uncontracted, unnormalized integrals of primitive Cartesian Gaussians
89 !> \param sphi_a contraction matrix for center a
90 !> \param sphi_b contraction matrix for center b
91 !> \param sphi_c contraction matrix for center c
92 !> \param ncoa number of cartesian orbitals on a
93 !> \param ncob number of cartesian orbitals on b
94 !> \param ncoc number of cartesian orbitals on c
95 !> \param nsgfa number of spherical Gaussian functions on a
96 !> \param nsgfb number of spherical Gaussian functions on b
97 !> \param nsgfc number of spherical Gaussian functions on c
98 ! **************************************************************************************************
99  SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
100  nsgfa, nsgfb, nsgfc)
101 
102  REAL(kind=dp), DIMENSION(:, :, :) :: abcint, sabc
103  REAL(kind=dp), DIMENSION(:, :) :: sphi_a, sphi_b, sphi_c
104  INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
105 
106  CHARACTER(LEN=*), PARAMETER :: routinen = 'abc_contract'
107 
108  INTEGER :: handle, i, m1, m2, m3, msphia, msphib, &
109  msphic
110  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpc, cpp
111 
112  CALL timeset(routinen, handle)
113 
114  cpassert(SIZE(abcint, 1) == nsgfa)
115  cpassert(SIZE(abcint, 2) == nsgfb)
116 
117  msphia = SIZE(sphi_a, 1)
118  msphib = SIZE(sphi_b, 1)
119  msphic = SIZE(sphi_c, 1)
120 
121  m1 = SIZE(sabc, 1)
122  m2 = SIZE(sabc, 2)
123  m3 = SIZE(sabc, 3)
124 
125  ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, nsgfb, m3))
126 
127  CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, cpp, nsgfa)
128  DO i = 1, m3
129  CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpp(:, :, i), nsgfa, sphi_b, msphib, &
130  0.0_dp, cpc(:, :, i), nsgfa)
131  END DO
132  CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, 1._dp, cpc, nsgfa*nsgfb, sphi_c, msphic, 0.0_dp, &
133  abcint, nsgfa*nsgfb)
134 
135  DEALLOCATE (cpp, cpc)
136 
137  CALL timestop(handle)
138 
139  END SUBROUTINE abc_contract
140 
141 ! **************************************************************************************************
142 !> \brief contract four-center overlap integrals (a,b,c,d) and transfer
143 !> to spherical Gaussians
144 !> \param abcdint contracted, normalized integrals of spherical Gaussians
145 !> \param sabcd uncontracted, unnormalized integrals of primitive Cartesian Gaussians
146 !> \param sphi_a contraction matrix for center a
147 !> \param sphi_b contraction matrix for center b
148 !> \param sphi_c contraction matrix for center c
149 !> \param sphi_d contraction matrix for center d
150 !> \param ncoa number of cartesian orbitals on a
151 !> \param ncob number of cartesian orbitals on b
152 !> \param ncoc number of cartesian orbitals on c
153 !> \param ncod number of cartesian orbitals on d
154 !> \param nsgfa number of spherical Gaussian functions on a
155 !> \param nsgfb number of spherical Gaussian functions on b
156 !> \param nsgfc number of spherical Gaussian functions on c
157 !> \param nsgfd number of spherical Gaussian functions on d
158 ! **************************************************************************************************
159  SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, &
160  ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
161 
162  REAL(kind=dp), DIMENSION(:, :, :, :), &
163  INTENT(INOUT) :: abcdint
164  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sabcd
165  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sphi_a, sphi_b, sphi_c, sphi_d
166  INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, &
167  nsgfc, nsgfd
168 
169  CHARACTER(LEN=*), PARAMETER :: routinen = 'abcd_contract'
170 
171  INTEGER :: handle, isgfc, isgfd, m1, m2, m3, m4, &
172  msphia, msphib, msphic, msphid
173  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: temp_cccc, work_cpcc
174  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: temp_cpcc, work_cppc
175  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: cpcc, cppc, cppp
176 
177  CALL timeset(routinen, handle)
178 
179  msphia = SIZE(sphi_a, 1)
180  msphib = SIZE(sphi_b, 1)
181  msphic = SIZE(sphi_c, 1)
182  msphid = SIZE(sphi_d, 1)
183 
184  m1 = SIZE(sabcd, 1)
185  m2 = SIZE(sabcd, 2)
186  m3 = SIZE(sabcd, 3)
187  m4 = SIZE(sabcd, 4)
188 
189  ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
190  cpcc(nsgfa, m2, nsgfc, nsgfd))
191 
192  ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
193  work_cppc = 0._dp
194  temp_cpcc = 0._dp
195 
196  ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
197  work_cpcc = 0._dp
198  temp_cccc = 0._dp
199 
200  CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
201  0.0_dp, cppp, nsgfa)
202  CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
203  sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
204 
205  DO isgfd = 1, nsgfd
206  work_cppc(:, :, :) = cppc(:, :, :, isgfd)
207  CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
208  sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
209  cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
210  DO isgfc = 1, nsgfc
211  work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
212  CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
213  msphib, 0.0_dp, temp_cccc, nsgfa)
214  abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
215  END DO
216  END DO
217 
218  DEALLOCATE (cpcc, cppc, cppp)
219  DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
220 
221  CALL timestop(handle)
222 
223  END SUBROUTINE abcd_contract
224 
225 ! **************************************************************************************************
226 !> \brief 3-center contraction routine from primitive cartesain Gaussians to spherical Gaussian
227 !> functions. Exploits LIBXSMM for performance, falls back to BLAS if LIBXSMM not available.
228 !> Requires pre-allocation of work buffers and pre-transposition of the sphi_a array. Requires
229 !> the LIBXSMM library to be initialized somewhere before this routine is called.
230 !> \param abcint contracted integrals
231 !> \param sabc uncontracted integrals
232 !> \param tsphi_a assumed to have dimensions nsgfa x ncoa
233 !> \param sphi_b assumed to have dimensions ncob x nsgfb
234 !> \param sphi_c assumed to have dimensions ncoc x nsgfc
235 !> \param ncoa ...
236 !> \param ncob ...
237 !> \param ncoc ...
238 !> \param nsgfa ...
239 !> \param nsgfb ...
240 !> \param nsgfc ...
241 !> \param cpp_buffer ...
242 !> \param ccp_buffer ...
243 !> \note tested from version 1.9.0 of libxsmm
244 ! **************************************************************************************************
245  SUBROUTINE libxsmm_abc_contract(abcint, sabc, tsphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
246  nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer)
247 
248  REAL(dp), DIMENSION(*) :: abcint, sabc
249  REAL(dp), DIMENSION(:, :) :: tsphi_a, sphi_b, sphi_c
250  INTEGER, INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
251  REAL(dp), DIMENSION(nsgfa*ncob) :: cpp_buffer
252  REAL(dp), DIMENSION(nsgfa*nsgfb*ncoc) :: ccp_buffer
253 
254  CHARACTER(LEN=*), PARAMETER :: routinen = 'libxsmm_abc_contract'
255 
256  INTEGER :: handle, i
257  LOGICAL :: libxsmm_kernels_available
258 #if(__LIBXSMM)
259  INTEGER(libxsmm_blasint_kind) :: m, n, k
260  TYPE(libxsmm_dmmfunction) :: xmm1, xmm2
261 #endif
262 
263  CALL timeset(routinen, handle)
264  libxsmm_kernels_available = .false.
265 
266 #if(__LIBXSMM)
267  !We make use of libxsmm code dispatching feature and call the same kernel multiple times
268  !We loop over the last index of the matrix and call libxsmm each time
269 
270  !pre-fetch the kernels
271  m = nsgfa; n = ncob; k = ncoa
272  CALL libxsmm_dispatch(xmm1, m, n, k, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
273  m = nsgfa; n = nsgfb; k = ncob
274  CALL libxsmm_dispatch(xmm2, m, n, k, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
275 
276  libxsmm_kernels_available = libxsmm_available(xmm1) .AND. libxsmm_available(xmm2)
277 
278  IF (libxsmm_kernels_available) THEN
279  ! contractions over a and b
280  DO i = 1, ncoc
281  CALL libxsmm_dmmcall(xmm1, tsphi_a, sabc((i - 1)*ncoa*ncob + 1), cpp_buffer)
282  CALL libxsmm_dmmcall(xmm2, cpp_buffer, sphi_b, ccp_buffer((i - 1)*nsgfa*nsgfb + 1))
283  END DO
284  ELSE
285  cpwarn("libxsmm available, but kernels are not, fallback to dgemm")
286  END IF
287 #endif
288 
289  IF (.NOT. libxsmm_kernels_available) THEN
290  ! we follow the same flow as for libxsmm above, but use BLAS dgemm
291  ! contractions over a and b
292  DO i = 1, ncoc
293  CALL dgemm("N", "N", nsgfa, ncob, ncoa, 1.0_dp, tsphi_a, nsgfa, sabc((i - 1)*ncoa*ncob + 1), &
294  ncoa, 0.0_dp, cpp_buffer, nsgfa)
295  CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1.0_dp, cpp_buffer, nsgfa, sphi_b, ncob, 0.0_dp, &
296  ccp_buffer((i - 1)*nsgfa*nsgfb + 1), nsgfa)
297  END DO
298  END IF
299 
300  ! last contraction, over c, as a larger MM
301  CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, 1.0_dp, ccp_buffer, nsgfa*nsgfb, sphi_c, ncoc, &
302  0.0_dp, abcint, nsgfa*nsgfb)
303 
304  CALL timestop(handle)
305 
306  END SUBROUTINE libxsmm_abc_contract
307 
308 END MODULE ai_contraction_sphi
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.
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
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
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
subroutine, public libxsmm_abc_contract(abcint, sabc, tsphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer)
3-center contraction routine from primitive cartesain Gaussians to spherical Gaussian functions....
subroutine, public ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
contract overlap integrals (a,b) and transfer to spherical Gaussians
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34