21 double *
const coef_xyz) {
23 assert(coef_xyz != NULL);
26 const int lp = (coef->
size[0] - 1);
27 for (
int lzp = 0; lzp <= lp; lzp++) {
28 for (
int lyp = 0; lyp <= lp - lzp; lyp++) {
29 for (
int lxp = 0; lxp <= lp - lzp - lyp; lxp++, lxyz++) {
30 coef_xyz[lxyz] =
idx3(coef[0], lzp, lyp, lxp);
37 double *
const coef_xyz) {
39 assert(coef_xyz != NULL);
41 const int lp = (coef->
size[0] - 1);
42 for (
int lzp = 0; lzp <= lp; lzp++) {
43 for (
int lyp = 0; lyp <= lp - lzp; lyp++) {
44 for (
int lxp = 0; lxp <= lp - lzp - lyp; lxp++, lxyz++) {
45 coef_xyz[lxyz] =
idx3(coef[0], lyp, lxp, lzp);
54 assert(coef_xyz != NULL);
56 const int lp = coef->
size[0] - 1;
57 for (
int lzp = 0; lzp <= lp; lzp++) {
58 for (
int lyp = 0; lyp <= lp - lzp; lyp++) {
59 for (
int lxp = 0; lxp <= lp - lzp - lyp; lxp++, lxyz++) {
60 idx3(coef[0], lzp, lyp, lxp) = coef_xyz[lxyz];
63 for (
int lxp = lp - lzp - lyp + 1; lxp <= lp; lxp++)
64 idx3(coef[0], lzp, lyp, lxp) = 0.0;
73 const int *lmin,
const int *lmax,
const int lp,
const double prefactor,
81 assert(alpha != NULL);
82 assert(coef_xyz != NULL);
83 assert(coef_xyz->
data != NULL);
87 for (
int lzb = 0; lzb <= lmax[1]; lzb++) {
88 for (
int lyb = 0; lyb <= lmax[1] - lzb; lyb++) {
89 const int lxb_min =
imax(lmin[1] - lzb - lyb, 0);
90 for (
int lxb = lxb_min; lxb <= lmax[1] - lzb - lyb; lxb++) {
91 const int jco = coset(lxb, lyb, lzb);
92 for (
int lza = 0; lza <= lmax[0]; lza++) {
93 for (
int lya = 0; lya <= lmax[0] - lza; lya++) {
94 const int lxa_min =
imax(lmin[0] - lza - lya, 0);
95 for (
int lxa = lxa_min; lxa <= lmax[0] - lza - lya; lxa++) {
96 const int ico = coset(lxa, lya, lza);
97 const double pab_ =
idx2(pab[0], jco, ico);
98 for (
int lxp = 0; lxp <= lxa + lxb; lxp++) {
100 idx4(alpha[0], 0, lxb, lxa, lxp) * pab_ * prefactor;
101 for (
int lzp = 0; lzp <= lp - lxa - lxb; lzp++) {
102 for (
int lyp = 0; lyp <= lp - lxa - lxb - lzp; lyp++) {
103 const double p2 =
idx4(alpha[0], 1, lyb, lya, lyp) *
104 idx4(alpha[0], 2, lzb, lza, lzp) * p1;
105 idx3(coef_xyz[0], lxp, lzp, lyp) += p2;
121 const int *
const lmin,
const int *
const lmax,
const int lp,
122 const double prefactor,
123 const tensor *
const alpha,
129 assert(alpha != NULL);
130 assert(coef_xyz != NULL);
131 assert(coef_xyz->
data != NULL);
136 for (
int lzb = 0; lzb <= lmax[1]; lzb++) {
137 for (
int lyb = 0; lyb <= lmax[1] - lzb; lyb++) {
138 const int lxb_min =
imax(lmin[1] - lzb - lyb, 0);
139 for (
int lxb = lxb_min; lxb <= lmax[1] - lzb - lyb; lxb++) {
140 const int jco = coset(lxb, lyb, lzb);
141 for (
int lza = 0; lza <= lmax[0]; lza++) {
142 for (
int lya = 0; lya <= lmax[0] - lza; lya++) {
143 const int lxa_min =
imax(lmin[0] - lza - lya, 0);
144 for (
int lxa = lxa_min; lxa <= lmax[0] - lza - lya; lxa++) {
145 const int ico = coset(lxa, lya, lza);
152 for (
int lxp = 0; lxp <= lxa + lxb; lxp++) {
153 for (
int lzp = 0; lzp <= lp - lxa - lxb; lzp++) {
154 for (
int lyp = 0; lyp <= lp - lxa - lxb - lzp; lyp++) {
155 const double p2 =
idx4(alpha[0], 1, lyb, lya, lyp) *
156 idx4(alpha[0], 2, lzb, lza, lzp) *
157 idx4(alpha[0], 0, lxb, lxa, lxp) *
159 pab_ +=
idx3(coef_xyz[0], lyp, lxp, lzp) * p2;
163 idx2(vab[0], jco, ico) += pab_;
174 const double rp[3],
const int *lmax,
176 assert(alpha != NULL);
178 memset(alpha->data, 0, alpha->alloc_size_ *
sizeof(
double));
185 for (
int iaxis = 0; iaxis < 3; iaxis++) {
186 const double drpa = rp[iaxis] - ra[iaxis];
187 const double drpb = rp[iaxis] - rb[iaxis];
188 for (
int lxa = 0; lxa <= lmax[0]; lxa++) {
189 for (
int lxb = 0; lxb <= lmax[1]; lxb++) {
190 double binomial_k_lxa = 1.0;
192 for (
int k = 0; k <= lxa; k++) {
193 double binomial_l_lxb = 1.0;
195 for (
int l = 0; l <= lxb; l++) {
196 idx4(alpha[0], iaxis, lxb, lxa, lxa - l + lxb - k) +=
197 binomial_k_lxa * binomial_l_lxb * a * b;
198 binomial_l_lxb *= ((double)(lxb - l)) / ((
double)(l + 1));
201 binomial_k_lxa *= ((double)(lxa - k)) / ((
double)(k + 1));
215 assert(coef_xyz != NULL);
216 const int lp = coef_xyz->
size[0] - 1;
231 if (coef_ijk.
data == NULL)
241 for (
int i = 0;
i < 3;
i++) {
242 for (
int j = 0; j < 3; j++) {
243 idx3(hmatgridp, 0, j,
i) = 1.0;
244 for (
int k = 1; k <= lp; k++) {
245 idx3(hmatgridp, k, j,
i) =
idx3(hmatgridp, k - 1, j,
i) * dh[j][
i];
251 for (
int klx = 0; klx <= lpx; klx++) {
252 for (
int jlx = 0; jlx <= lpx - klx; jlx++) {
253 for (
int ilx = 0; ilx <= lpx - klx - jlx; ilx++) {
254 const int lx = ilx + jlx + klx;
255 const int lpy = lp - lx;
256 const double tx =
idx3(hmatgridp, ilx, 0, 0) *
257 idx3(hmatgridp, jlx, 1, 0) *
258 idx3(hmatgridp, klx, 2, 0) * fac(lx) *
inv_fac[klx] *
261 for (
int kly = 0; kly <= lpy; kly++) {
262 for (
int jly = 0; jly <= lpy - kly; jly++) {
263 for (
int ily = 0; ily <= lpy - kly - jly; ily++) {
264 const int ly = ily + jly + kly;
265 const int lpz = lp - lx - ly;
266 const double ty = tx *
idx3(hmatgridp, ily, 0, 1) *
267 idx3(hmatgridp, jly, 1, 1) *
268 idx3(hmatgridp, kly, 2, 1) * fac(ly) *
270 for (
int klz = 0; klz <= lpz; klz++) {
271 for (
int jlz = 0; jlz <= lpz - klz; jlz++) {
272 for (
int ilz = 0; ilz <= lpz - klz - jlz; ilz++) {
273 const int lz = ilz + jlz + klz;
274 const int il = ilx + ily + ilz;
275 const int jl = jlx + jly + jlz;
276 const int kl = klx + kly + klz;
281 idx3(coef_ijk, il, kl, jl) +=
282 idx3(coef_xyz[0], lx, lz, ly) * ty *
283 idx3(hmatgridp, ilz, 0, 2) *
284 idx3(hmatgridp, jlz, 1, 2) *
285 idx3(hmatgridp, klz, 2, 2) * fac(lz) *
inv_fac[klz] *
308 const int lp = coef_xyz->
size[0] - 1;
322 if (coef_ijk.
data == NULL)
332 for (
int i = 0;
i < 3;
i++) {
333 for (
int j = 0; j < 3; j++) {
334 idx3(hmatgridp, 0, j,
i) = 1.0;
335 for (
int k = 1; k <= lp; k++) {
336 idx3(hmatgridp, k, j,
i) =
idx3(hmatgridp, k - 1, j,
i) * dh[j][
i];
342 for (
int klx = 0; klx <= lpx; klx++) {
343 for (
int jlx = 0; jlx <= lpx - klx; jlx++) {
344 for (
int ilx = 0; ilx <= lpx - klx - jlx; ilx++) {
345 const int lx = ilx + jlx + klx;
346 const int lpy = lp - lx;
347 for (
int kly = 0; kly <= lpy; kly++) {
348 for (
int jly = 0; jly <= lpy - kly; jly++) {
349 for (
int ily = 0; ily <= lpy - kly - jly; ily++) {
350 const int ly = ily + jly + kly;
351 const int lpz = lp - lx - ly;
352 for (
int klz = 0; klz <= lpz; klz++) {
353 for (
int jlz = 0; jlz <= lpz - klz; jlz++) {
354 for (
int ilz = 0; ilz <= lpz - klz - jlz; ilz++) {
355 const int lz = ilz + jlz + klz;
356 const int il = ilx + ily + ilz;
357 const int jl = jlx + jly + jlz;
358 const int kl = klx + kly + klz;
363 idx3(coef_ijk, ly, lx, lz) +=
364 idx3(coef_xyz[0], jl, il, kl) *
365 idx3(hmatgridp, ilx, 0, 0) *
366 idx3(hmatgridp, jlx, 1, 0) *
367 idx3(hmatgridp, klx, 2, 0) *
368 idx3(hmatgridp, ily, 0, 1) *
369 idx3(hmatgridp, jly, 1, 1) *
370 idx3(hmatgridp, kly, 2, 1) *
371 idx3(hmatgridp, ilz, 0, 2) *
372 idx3(hmatgridp, jlz, 1, 2) *
373 idx3(hmatgridp, klz, 2, 2) * fac(lx) * fac(ly) *
375 (fac(ilx) * fac(ily) * fac(ilz) * fac(jlx) * fac(jly) *
376 fac(jlz) * fac(klx) * fac(kly) * fac(klz));
void grid_compute_coefficients_dgemm(const int *lmin, const int *lmax, const int lp, const double prefactor, const tensor *const alpha, const tensor *const pab, tensor *coef_xyz)
void grid_compute_vab(const int *const lmin, const int *const lmax, const int lp, const double prefactor, const tensor *const alpha, const tensor *const coef_xyz, tensor *vab)