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