(git:b279b6b)
grid_hip_prepare_pab.h
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: BSD-3-Clause */
6 /*----------------------------------------------------------------------------*/
7 
8 /*
9  * NOTE : derived from the reference and GPU backends
10  * Authors :
11  - Dr Mathieu Taillefumier (ETH Zurich / CSCS)
12  - Advanced Micro Devices, Inc.
13 */
14 
15 #include "../common/grid_constants.h"
16 #include <stdbool.h>
17 #include <stdio.h>
18 #include <stdlib.h>
19 
20 namespace rocm_backend {
21 /*******************************************************************************
22  * \brief Implementation of function GRID_FUNC_AB, ie. identity transformation.
23  ******************************************************************************/
24 template <typename T>
25 __device__ __inline__ void prepare_pab_AB(const orbital a, const orbital b,
26  const T pab_val, const int n,
27  T *cab) {
28 
29  // simply copy pab to cab
30  prep_term<T>(a, b, pab_val, n, cab);
31 }
32 
33 /*******************************************************************************
34  * \brief Implementation of function GRID_FUNC_DADB.
35  ******************************************************************************/
36 template <typename T>
37 __device__ __inline__ void
38 prepare_pab_DADB(const orbital a, const orbital b, const T zeta, const T zetb,
39  const T pab_val, const int n, T *cab) {
40 
41  // creates cab such that mapping it with pgf_a pgf_b
42  // is equivalent to mapping pab with 0.5 * (nabla pgf_a) . (nabla pgf_b)
43  // (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx
44  // pgf_{b-1x} - 2*zetb*pgf_{b+1x})
45 
46  for (int i = 0; i < 3; i++) {
47  prep_term<T>(down(i, a), down(i, b), 0.5 * a.l[i] * b.l[i] * pab_val, n,
48  cab);
49  prep_term<T>(down(i, a), up(i, b), -1.0 * a.l[i] * zetb * pab_val, n, cab);
50  prep_term<T>(up(i, a), down(i, b), -1.0 * zeta * b.l[i] * pab_val, n, cab);
51  prep_term<T>(up(i, a), up(i, b), 2.0 * zeta * zetb * pab_val, n, cab);
52  }
53 }
54 
55 /*******************************************************************************
56  * \brief Implementation of function GRID_FUNC_ADBmDAB_{X,Y,Z}.
57  ******************************************************************************/
58 template <typename T>
59 __device__ void prepare_pab_ADBmDAB(const int idir, const orbital a,
60  const orbital b, const T zeta, const T zetb,
61  const T pab_val, const int n, T *cab) {
62 
63  // creates cab such that mapping it with pgf_a pgf_b
64  // is equivalent to mapping pab with
65  // pgf_a (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) pgf_b
66  // ( pgf_a ) (ddx pgf_b) - (ddx pgf_a)( pgf_b ) =
67  // pgf_a *(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) -
68  // (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b
69 
70  prep_term<T>(a, down(idir, b), +b.l[idir] * pab_val, n, cab);
71  prep_term<T>(a, up(idir, b), -2.0 * zetb * pab_val, n, cab);
72  prep_term<T>(down(idir, a), b, -a.l[idir] * pab_val, n, cab);
73  prep_term<T>(up(idir, a), b, +2.0 * zeta * pab_val, n, cab);
74 }
75 
76 /*******************************************************************************
77  * \brief Implementation of function GRID_FUNC_ARDBmDARB_{X,Y,Z}{X,Y,Z}.
78  ******************************************************************************/
79 template <typename T>
80 __device__ void prepare_pab_ARDBmDARB(const int idir, const int ir,
81  const orbital a, const orbital b,
82  const T zeta, const T zetb,
83  const T pab_val, const int n, T *cab) {
84 
85  // creates cab such that mapping it with pgf_a pgf_b
86  // is equivalent to mapping pab with
87  // pgf_a (r-Rb)_{ir} (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) (r-Rb)_{ir}
88  // pgf_b ( pgf_a )(r-Rb)_{ir} (ddx pgf_b) - (ddx pgf_a) (r-Rb)_{ir} ( pgf_b )
89  // =
90  // pgf_a *(lbx pgf_{b-1x+1ir} - 2*zetb*pgf_{b+1x+1ir})
91  // -
92  // (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_{b+1ir}
93 
94  prep_term<T>(a, down(idir, up(ir, b)), b.l[idir] * pab_val, n, cab);
95  prep_term<T>(a, up(idir, up(ir, b)), -2.0 * zetb * pab_val, n, cab);
96  prep_term<T>(down(idir, a), up(ir, b), -a.l[idir] * pab_val, n, cab);
97  prep_term<T>(up(idir, a), up(ir, b), +2.0 * zeta * pab_val, n, cab);
98 }
99 
100 /*******************************************************************************
101  * \brief Implementation of function GRID_FUNC_DABpADB_{X,Y,Z}.
102  * \author Ole Schuett
103  ******************************************************************************/
104 template <typename T>
105 __device__ void prepare_pab_DABpADB(const int idir, const orbital a,
106  const orbital b, const T zeta, const T zetb,
107  const T pab_val, const int n, T *cab) {
108 
109  // creates cab such that mapping it with pgf_a pgf_b
110  // is equivalent to mapping pab with
111  // pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b
112  // ( pgf_a ) (ddx pgf_b) + (ddx pgf_a)( pgf_b ) =
113  // pgf_a *(lbx pgf_{b-1x} + 2*zetb*pgf_{b+1x}) +
114  // (lax pgf_{a-1x} + 2*zeta*pgf_{a+1x}) pgf_b
115 
116  prep_term<T>(a, down(idir, b), b.l[idir] * pab_val, n, cab);
117  prep_term<T>(a, up(idir, b), -2.0 * zetb * pab_val, n, cab);
118  prep_term<T>(down(idir, a), b, a.l[idir] * pab_val, n, cab);
119  prep_term<T>(up(idir, a), b, -2.0 * zeta * pab_val, n, cab);
120 }
121 
122 /*******************************************************************************
123  * \brief Implementation of function GRID_FUNC_{DX,DY,DZ}.
124  * \author Ole Schuett
125  ******************************************************************************/
126 template <typename T>
127 __device__ __inline__ void
128 prepare_pab_Di(const int ider, const orbital a, const orbital b, const T zeta,
129  const T zetb, const T pab_val, const int n, T *cab) {
130 
131  // creates cab such that mapping it with pgf_a pgf_b
132  // is equivalent to mapping pab with
133  // d_{ider1} pgf_a d_{ider1} pgf_b
134  // dx pgf_a dx pgf_b =
135  // (lax pgf_{a-1x})*(lbx pgf_{b-1x}) - 2*zetb*lax*pgf_{a-1x}*pgf_{b+1x}
136  // -
137  // lbx pgf_{b-1x}*2*zeta*pgf_{a+1x}+ 4*zeta*zetab*pgf_{a+1x}pgf_{b+1x}
138 
139  prep_term<T>(down(ider, a), down(ider, b), a.l[ider] * b.l[ider] * pab_val, n,
140  cab);
141  prep_term<T>(down(ider, a), up(ider, b), -2.0 * a.l[ider] * zetb * pab_val, n,
142  cab);
143  prep_term<T>(up(ider, a), down(ider, b), -2.0 * zeta * b.l[ider] * pab_val, n,
144  cab);
145  prep_term<T>(up(ider, a), up(ider, b), +4.0 * zeta * zetb * pab_val, n, cab);
146 }
147 
148 /*******************************************************************************
149  * \brief Helper for grid_prepare_pab_DiDj.
150  ******************************************************************************/
151 template <typename T>
152 __device__ __inline__ void oneterm_dijdij(const int idir, const T func_a,
153  const orbital a, const orbital b,
154  const T zetb, const int n, T *cab) {
155 
156  int i1, i2;
157  if (idir == 0) {
158  i1 = 0;
159  i2 = 1;
160  } else if (idir == 1) {
161  i1 = 1;
162  i2 = 2;
163  } else if (idir == 2) {
164  i1 = 2;
165  i2 = 0;
166  } else {
167  return; // error
168  }
169 
170  prep_term<T>(a, down(i1, down(i2, b)), b.l[i1] * b.l[i2] * func_a, n, cab);
171  prep_term<T>(a, up(i1, down(i2, b)), -2.0 * zetb * b.l[i2] * func_a, n, cab);
172  prep_term<T>(a, down(i1, up(i2, b)), -2.0 * zetb * b.l[i1] * func_a, n, cab);
173  prep_term<T>(a, up(i1, up(i2, b)), +4.0 * zetb * zetb * func_a, n, cab);
174 }
175 
176 /*******************************************************************************
177  * \brief Implementation of function GRID_FUNC_{DXDY,DYDZ,DZDX}
178  * \author Ole Schuett
179  ******************************************************************************/
180 template <typename T>
181 __device__ __inline__ void
182 prepare_pab_DiDj(const int ider1, const int ider2, const orbital a,
183  const orbital b, const T zeta, const T zetb, const T pab_val,
184  const int n, T *cab) {
185 
186  // creates cab such that mapping it with pgf_a pgf_b
187  // is equivalent to mapping pab with
188  // d_{ider1} pgf_a d_{ider1} pgf_b
189 
190  const T func_a1 = a.l[ider1] * a.l[ider2] * pab_val;
191  oneterm_dijdij<T>(ider1, func_a1, down(ider1, down(ider2, a)), b, zetb, n,
192  cab);
193 
194  const T func_a2 = -2.0 * zeta * a.l[ider2] * pab_val;
195  oneterm_dijdij<T>(ider1, func_a2, up(ider1, down(ider2, a)), b, zetb, n, cab);
196 
197  const T func_a3 = -2.0 * zeta * a.l[ider1] * pab_val;
198  oneterm_dijdij<T>(ider1, func_a3, down(ider1, up(ider2, a)), b, zetb, n, cab);
199 
200  const T func_a4 = 4.0 * zeta * zeta * pab_val;
201  oneterm_dijdij<T>(ider1, func_a4, up(ider1, up(ider2, a)), b, zetb, n, cab);
202 }
203 
204 /*******************************************************************************
205  * \brief Helper for grid_prepare_pab_Di2.
206  * \author Ole Schuett
207  ******************************************************************************/
208 template <typename T>
209 __device__ void oneterm_diidii(const int idir, const T func_a, const orbital a,
210  const orbital b, const T zetb, const int n,
211  T *cab) {
212 
213  prep_term<T>(a, down(idir, down(idir, b)),
214  b.l[idir] * (b.l[idir] - 1) * func_a, n, cab);
215  prep_term<T>(a, b, -2.0 * zetb * (2 * b.l[idir] + 1) * func_a, n, cab);
216  prep_term<T>(a, up(idir, up(idir, b)), +4.0 * zetb * zetb * func_a, n, cab);
217 }
218 
219 /*******************************************************************************
220  * \brief Implementation of function GRID_FUNC_{DXDX,DYDY,DZDZ}
221  ******************************************************************************/
222 template <typename T>
223 __device__ __inline__ void
224 prepare_pab_Di2(const int ider, const orbital a, const orbital b, const T zeta,
225  const T zetb, const T pab_val, const int n, T *cab) {
226 
227  // creates cab such that mapping it with pgf_a pgf_b
228  // is equivalent to mapping pab with
229  // dd_{ider1} pgf_a dd_{ider1} pgf_b
230 
231  const T func_a1 = a.l[ider] * (a.l[ider] - 1) * pab_val;
232  oneterm_diidii<T>(ider, func_a1, down(ider, down(ider, a)), b, zetb, n, cab);
233 
234  const T func_a2 = -2.0 * zeta * (2 * a.l[ider] + 1) * pab_val;
235  oneterm_diidii<T>(ider, func_a2, a, b, zetb, n, cab);
236 
237  const T func_a3 = 4.0 * zeta * zeta * pab_val;
238  oneterm_diidii<T>(ider, func_a3, up(ider, up(ider, a)), b, zetb, n, cab);
239 }
240 
241 /*******************************************************************************
242  * \brief Transforms a given element of the density matrix according to func.
243  * \param func Transformation function to apply, one of GRID_FUNC_*.
244  * \param {a,b} Orbital angular momenta.
245  * \param zet_{a,b} Gaussian exponents.
246  * \param pab_val Input matrix element of pab.
247  * \param n Leading dimensions of output matrix cab.
248  * \param cab Output matrix.
249  ******************************************************************************/
250 template <typename T>
251 __device__ __inline__ void
252 prepare_pab(const enum grid_func func, const orbital a, const orbital b,
253  const T zeta, const T zetb, const T pab_val, const int n, T *cab) {
254 
255  // This switch statment will be in an inner loop but only with few iterations.
256  switch (func) {
257  case GRID_FUNC_AB:
258  prepare_pab_AB<T>(a, b, pab_val, n, cab);
259  break;
260  case GRID_FUNC_DADB:
261  prepare_pab_DADB<T>(a, b, zeta, zetb, pab_val, n, cab);
262  break;
263  case GRID_FUNC_ADBmDAB_X:
264  prepare_pab_ADBmDAB<T>(0, a, b, zeta, zetb, pab_val, n, cab);
265  break;
266  case GRID_FUNC_ADBmDAB_Y:
267  prepare_pab_ADBmDAB<T>(1, a, b, zeta, zetb, pab_val, n, cab);
268  break;
269  case GRID_FUNC_ADBmDAB_Z:
270  prepare_pab_ADBmDAB<T>(2, a, b, zeta, zetb, pab_val, n, cab);
271  break;
273  prepare_pab_ARDBmDARB<T>(0, 0, a, b, zeta, zetb, pab_val, n, cab);
274  break;
276  prepare_pab_ARDBmDARB<T>(0, 1, a, b, zeta, zetb, pab_val, n, cab);
277  break;
279  prepare_pab_ARDBmDARB<T>(0, 2, a, b, zeta, zetb, pab_val, n, cab);
280  break;
282  prepare_pab_ARDBmDARB<T>(1, 0, a, b, zeta, zetb, pab_val, n, cab);
283  break;
285  prepare_pab_ARDBmDARB<T>(1, 1, a, b, zeta, zetb, pab_val, n, cab);
286  break;
288  prepare_pab_ARDBmDARB<T>(1, 2, a, b, zeta, zetb, pab_val, n, cab);
289  break;
291  prepare_pab_ARDBmDARB<T>(2, 0, a, b, zeta, zetb, pab_val, n, cab);
292  break;
294  prepare_pab_ARDBmDARB<T>(2, 1, a, b, zeta, zetb, pab_val, n, cab);
295  break;
297  prepare_pab_ARDBmDARB<T>(2, 2, a, b, zeta, zetb, pab_val, n, cab);
298  break;
299  case GRID_FUNC_DABpADB_X:
300  prepare_pab_DABpADB<T>(0, a, b, zeta, zetb, pab_val, n, cab);
301  break;
302  case GRID_FUNC_DABpADB_Y:
303  prepare_pab_DABpADB<T>(1, a, b, zeta, zetb, pab_val, n, cab);
304  break;
305  case GRID_FUNC_DABpADB_Z:
306  prepare_pab_DABpADB<T>(2, a, b, zeta, zetb, pab_val, n, cab);
307  break;
308  case GRID_FUNC_DX:
309  prepare_pab_Di<T>(0, a, b, zeta, zetb, pab_val, n, cab);
310  break;
311  case GRID_FUNC_DY:
312  prepare_pab_Di<T>(1, a, b, zeta, zetb, pab_val, n, cab);
313  break;
314  case GRID_FUNC_DZ:
315  prepare_pab_Di<T>(2, a, b, zeta, zetb, pab_val, n, cab);
316  break;
317  case GRID_FUNC_DXDY:
318  prepare_pab_DiDj<T>(0, 1, a, b, zeta, zetb, pab_val, n, cab);
319  break;
320  case GRID_FUNC_DYDZ:
321  prepare_pab_DiDj<T>(1, 2, a, b, zeta, zetb, pab_val, n, cab);
322  break;
323  case GRID_FUNC_DZDX:
324  prepare_pab_DiDj<T>(2, 0, a, b, zeta, zetb, pab_val, n, cab);
325  break;
326  case GRID_FUNC_DXDX:
327  prepare_pab_Di2<T>(0, a, b, zeta, zetb, pab_val, n, cab);
328  break;
329  case GRID_FUNC_DYDY:
330  prepare_pab_Di2<T>(1, a, b, zeta, zetb, pab_val, n, cab);
331  break;
332  case GRID_FUNC_DZDZ:
333  prepare_pab_Di2<T>(2, a, b, zeta, zetb, pab_val, n, cab);
334  break;
335  default:
336  break; // Error: Unknown ga_gb_function - do nothing.
337  }
338 }
339 
340 /*******************************************************************************
341  * \brief Returns difference in angular momentum range for given func.
342  * \author Ole Schuett
343  ******************************************************************************/
344 inline ldiffs_value prepare_get_ldiffs(const enum grid_func func) {
345  ldiffs_value ldiffs;
346 
347  switch (func) {
348  case GRID_FUNC_AB:
349  ldiffs.la_max_diff = 0;
350  ldiffs.la_min_diff = 0;
351  ldiffs.lb_max_diff = 0;
352  ldiffs.lb_min_diff = 0;
353  break;
354  case GRID_FUNC_DADB:
355  case GRID_FUNC_ADBmDAB_X:
356  case GRID_FUNC_ADBmDAB_Y:
357  case GRID_FUNC_ADBmDAB_Z:
358  case GRID_FUNC_DABpADB_X:
359  case GRID_FUNC_DABpADB_Y:
360  case GRID_FUNC_DABpADB_Z:
361  ldiffs.la_max_diff = +1;
362  ldiffs.la_min_diff = -1;
363  ldiffs.lb_max_diff = +1;
364  ldiffs.lb_min_diff = -1;
365  break;
375  ldiffs.la_max_diff = +1;
376  ldiffs.la_min_diff = -1;
377  ldiffs.lb_max_diff = +2; // this is legit
378  ldiffs.lb_min_diff = -1;
379  break;
380  case GRID_FUNC_DX:
381  case GRID_FUNC_DY:
382  case GRID_FUNC_DZ:
383  ldiffs.la_max_diff = +1;
384  ldiffs.la_min_diff = -1;
385  ldiffs.lb_max_diff = +1;
386  ldiffs.lb_min_diff = -1;
387  break;
388  case GRID_FUNC_DXDY:
389  case GRID_FUNC_DYDZ:
390  case GRID_FUNC_DZDX:
391  case GRID_FUNC_DXDX:
392  case GRID_FUNC_DYDY:
393  case GRID_FUNC_DZDZ:
394  ldiffs.la_max_diff = +2;
395  ldiffs.la_min_diff = -2;
396  ldiffs.lb_max_diff = +2;
397  ldiffs.lb_min_diff = -2;
398  break;
399  default:
400  fprintf(stderr, "Error: Unknown ga_gb_function %i.\n", func);
401  abort();
402  }
403 
404  return ldiffs;
405 }
406 } // namespace rocm_backend
grid_func
@ GRID_FUNC_ADBmDAB_X
@ GRID_FUNC_DZ
@ GRID_FUNC_DABpADB_X
@ GRID_FUNC_DY
@ GRID_FUNC_ARDBmDARB_YY
@ GRID_FUNC_DX
@ GRID_FUNC_ARDBmDARB_XY
@ GRID_FUNC_ARDBmDARB_ZY
@ GRID_FUNC_DYDY
@ GRID_FUNC_ARDBmDARB_XZ
@ GRID_FUNC_DZDX
@ GRID_FUNC_ARDBmDARB_YZ
@ GRID_FUNC_AB
@ GRID_FUNC_DYDZ
@ GRID_FUNC_ADBmDAB_Z
@ GRID_FUNC_DXDY
@ GRID_FUNC_DXDX
@ GRID_FUNC_ADBmDAB_Y
@ GRID_FUNC_DABpADB_Z
@ GRID_FUNC_ARDBmDARB_ZX
@ GRID_FUNC_DZDZ
@ GRID_FUNC_ARDBmDARB_YX
@ GRID_FUNC_ARDBmDARB_ZZ
@ GRID_FUNC_ARDBmDARB_XX
@ GRID_FUNC_DABpADB_Y
@ GRID_FUNC_DADB
static void const int const int i
real(dp), dimension(3) a
Definition: ai_eri_debug.F:31
real(dp), dimension(3) b
Definition: ai_eri_debug.F:31
__device__ void prepare_pab_ADBmDAB(const int idir, const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_ADBmDAB_{X,Y,Z}.
__device__ __inline__ void oneterm_dijdij(const int idir, const T func_a, const orbital a, const orbital b, const T zetb, const int n, T *cab)
Helper for grid_prepare_pab_DiDj.
__device__ __inline__ void prepare_pab_Di2(const int ider, const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_{DXDX,DYDY,DZDZ}.
__device__ __inline__ void prepare_pab_DiDj(const int ider1, const int ider2, const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_{DXDY,DYDZ,DZDX}.
__device__ __inline__ void prepare_pab_AB(const orbital a, const orbital b, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_AB, ie. identity transformation.
ldiffs_value prepare_get_ldiffs(const enum grid_func func)
Returns difference in angular momentum range for given func.
__device__ void prepare_pab_DABpADB(const int idir, const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_DABpADB_{X,Y,Z}.
__device__ __inline__ void prepare_pab_Di(const int ider, const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_{DX,DY,DZ}.
__device__ void prepare_pab_ARDBmDARB(const int idir, const int ir, const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_ARDBmDARB_{X,Y,Z}{X,Y,Z}.
__inline__ __device__ orbital down(const int i, const orbital &a)
Decrease i'th component of given orbital angular momentum.
__device__ __inline__ orbital up(const int i, const orbital &a)
Increase i'th component of given orbital angular momentum.
__device__ void oneterm_diidii(const int idir, const T func_a, const orbital a, const orbital b, const T zetb, const int n, T *cab)
Helper for grid_prepare_pab_Di2.
__device__ __inline__ void prepare_pab(const enum grid_func func, const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Transforms a given element of the density matrix according to func.
__device__ __inline__ void prepare_pab_DADB(const orbital a, const orbital b, const T zeta, const T zetb, const T pab_val, const int n, T *cab)
Implementation of function GRID_FUNC_DADB.
Differences in angular momentum.
Orbital angular momentum.