1 /*
2 * Copyright (C) 2013- Qiming Sun <osirpt.sun@gmail.com>
3 *
4 */
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <math.h>
9 #include <assert.h>
10 #include "config.h"
11 #include "cint_bas.h"
12 #include "rys_roots.h"
13 #include "misc.h"
14 #include "g2e.h"
15
16 #define DEF_GXYZ(type, G, GX, GY, GZ) \
17 type *GX = G; \
18 type *GY = G + envs->g_size; \
19 type *GZ = G + envs->g_size * 2
20
CINTinit_int2e_EnvVars(CINTEnvVars * envs,FINT * ng,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env)21 void CINTinit_int2e_EnvVars(CINTEnvVars *envs, FINT *ng, FINT *shls,
22 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env)
23 {
24 envs->natm = natm;
25 envs->nbas = nbas;
26 envs->atm = atm;
27 envs->bas = bas;
28 envs->env = env;
29 envs->shls = shls;
30
31 const FINT i_sh = shls[0];
32 const FINT j_sh = shls[1];
33 const FINT k_sh = shls[2];
34 const FINT l_sh = shls[3];
35 envs->i_l = bas(ANG_OF, i_sh);
36 envs->j_l = bas(ANG_OF, j_sh);
37 envs->k_l = bas(ANG_OF, k_sh);
38 envs->l_l = bas(ANG_OF, l_sh);
39 envs->x_ctr[0] = bas(NCTR_OF, i_sh);
40 envs->x_ctr[1] = bas(NCTR_OF, j_sh);
41 envs->x_ctr[2] = bas(NCTR_OF, k_sh);
42 envs->x_ctr[3] = bas(NCTR_OF, l_sh);
43 envs->nfi = (envs->i_l+1)*(envs->i_l+2)/2;
44 envs->nfj = (envs->j_l+1)*(envs->j_l+2)/2;
45 envs->nfk = (envs->k_l+1)*(envs->k_l+2)/2;
46 envs->nfl = (envs->l_l+1)*(envs->l_l+2)/2;
47 envs->nf = envs->nfi * envs->nfk * envs->nfl * envs->nfj;
48
49 envs->ri = env + atm(PTR_COORD, bas(ATOM_OF, i_sh));
50 envs->rj = env + atm(PTR_COORD, bas(ATOM_OF, j_sh));
51 envs->rk = env + atm(PTR_COORD, bas(ATOM_OF, k_sh));
52 envs->rl = env + atm(PTR_COORD, bas(ATOM_OF, l_sh));
53
54 envs->common_factor = (M_PI*M_PI*M_PI)*2/SQRTPI
55 * CINTcommon_fac_sp(envs->i_l) * CINTcommon_fac_sp(envs->j_l)
56 * CINTcommon_fac_sp(envs->k_l) * CINTcommon_fac_sp(envs->l_l);
57 if (env[PTR_EXPCUTOFF] == 0) {
58 envs->expcutoff = EXPCUTOFF;
59 } else {
60 // +1 to ensure accuracy. See comments in function CINT2e_loop_nopt
61 envs->expcutoff = MAX(MIN_EXPCUTOFF, env[PTR_EXPCUTOFF]) + 1;
62 }
63
64 envs->gbits = ng[GSHIFT];
65 envs->ncomp_e1 = ng[POS_E1];
66 envs->ncomp_e2 = ng[POS_E2];
67 envs->ncomp_tensor = ng[TENSOR];
68
69 envs->li_ceil = envs->i_l + ng[IINC];
70 envs->lj_ceil = envs->j_l + ng[JINC];
71 envs->lk_ceil = envs->k_l + ng[KINC];
72 envs->ll_ceil = envs->l_l + ng[LINC];
73 envs->nrys_roots =(envs->li_ceil + envs->lj_ceil
74 + envs->lk_ceil + envs->ll_ceil)/2 + 1;
75
76 assert(i_sh < SHLS_MAX);
77 assert(j_sh < SHLS_MAX);
78 assert(k_sh < SHLS_MAX);
79 assert(l_sh < SHLS_MAX);
80 assert(envs->i_l < ANG_MAX);
81 assert(envs->j_l < ANG_MAX);
82 assert(envs->k_l < ANG_MAX);
83 assert(envs->l_l < ANG_MAX);
84 assert(bas(ATOM_OF,i_sh) >= 0);
85 assert(bas(ATOM_OF,j_sh) >= 0);
86 assert(bas(ATOM_OF,k_sh) >= 0);
87 assert(bas(ATOM_OF,l_sh) >= 0);
88 assert(bas(ATOM_OF,i_sh) < natm);
89 assert(bas(ATOM_OF,j_sh) < natm);
90 assert(bas(ATOM_OF,k_sh) < natm);
91 assert(bas(ATOM_OF,l_sh) < natm);
92 assert(envs->nrys_roots < MXRYSROOTS);
93
94 FINT dli, dlj, dlk, dll;
95 FINT ibase = envs->li_ceil > envs->lj_ceil;
96 FINT kbase = envs->lk_ceil > envs->ll_ceil;
97 if (envs->nrys_roots <= 2) { // use the fully optimized lj_4d algorithm
98 ibase = 0;
99 kbase = 0;
100 }
101 if (kbase) {
102 dlk = envs->lk_ceil + envs->ll_ceil + 1;
103 dll = envs->ll_ceil + 1;
104 } else {
105 dlk = envs->lk_ceil + 1;
106 dll = envs->lk_ceil + envs->ll_ceil + 1;
107 }
108
109 if (ibase) {
110 dli = envs->li_ceil + envs->lj_ceil + 1;
111 dlj = envs->lj_ceil + 1;
112 } else {
113 dli = envs->li_ceil + 1;
114 dlj = envs->li_ceil + envs->lj_ceil + 1;
115 }
116 envs->g_stride_i = envs->nrys_roots;
117 envs->g_stride_k = envs->nrys_roots * dli;
118 envs->g_stride_l = envs->nrys_roots * dli * dlk;
119 envs->g_stride_j = envs->nrys_roots * dli * dlk * dll;
120 envs->g_size = envs->nrys_roots * dli * dlk * dll * dlj;
121
122 if (kbase) {
123 envs->g2d_klmax = envs->g_stride_k;
124 envs->rx_in_rklrx = envs->rk;
125 envs->rkrl[0] = envs->rk[0] - envs->rl[0];
126 envs->rkrl[1] = envs->rk[1] - envs->rl[1];
127 envs->rkrl[2] = envs->rk[2] - envs->rl[2];
128 } else {
129 envs->g2d_klmax = envs->g_stride_l;
130 envs->rx_in_rklrx = envs->rl;
131 envs->rkrl[0] = envs->rl[0] - envs->rk[0];
132 envs->rkrl[1] = envs->rl[1] - envs->rk[1];
133 envs->rkrl[2] = envs->rl[2] - envs->rk[2];
134 }
135
136 if (ibase) {
137 envs->g2d_ijmax = envs->g_stride_i;
138 envs->rx_in_rijrx = envs->ri;
139 envs->rirj[0] = envs->ri[0] - envs->rj[0];
140 envs->rirj[1] = envs->ri[1] - envs->rj[1];
141 envs->rirj[2] = envs->ri[2] - envs->rj[2];
142 } else {
143 envs->g2d_ijmax = envs->g_stride_j;
144 envs->rx_in_rijrx = envs->rj;
145 envs->rirj[0] = envs->rj[0] - envs->ri[0];
146 envs->rirj[1] = envs->rj[1] - envs->ri[1];
147 envs->rirj[2] = envs->rj[2] - envs->ri[2];
148 }
149
150 if (kbase) {
151 if (ibase) {
152 envs->f_g0_2d4d = &CINTg0_2e_ik2d4d;
153 } else {
154 envs->f_g0_2d4d = &CINTg0_2e_kj2d4d;
155 }
156 } else {
157 if (ibase) {
158 envs->f_g0_2d4d = &CINTg0_2e_il2d4d;
159 } else {
160 envs->f_g0_2d4d = &CINTg0_2e_lj2d4d;
161 }
162 }
163 envs->f_g0_2e = &CINTg0_2e;
164 }
165
CINTg2e_index_xyz(FINT * idx,const CINTEnvVars * envs)166 void CINTg2e_index_xyz(FINT *idx, const CINTEnvVars *envs)
167 {
168 const FINT i_l = envs->i_l;
169 const FINT j_l = envs->j_l;
170 const FINT k_l = envs->k_l;
171 const FINT l_l = envs->l_l;
172 const FINT nfi = envs->nfi;
173 const FINT nfj = envs->nfj;
174 const FINT nfk = envs->nfk;
175 const FINT nfl = envs->nfl;
176 const FINT di = envs->g_stride_i;
177 const FINT dk = envs->g_stride_k;
178 const FINT dl = envs->g_stride_l;
179 const FINT dj = envs->g_stride_j;
180 FINT i, j, k, l, n;
181 FINT ofx, ofkx, oflx;
182 FINT ofy, ofky, ofly;
183 FINT ofz, ofkz, oflz;
184 FINT i_nx[CART_MAX], i_ny[CART_MAX], i_nz[CART_MAX];
185 FINT j_nx[CART_MAX], j_ny[CART_MAX], j_nz[CART_MAX];
186 FINT k_nx[CART_MAX], k_ny[CART_MAX], k_nz[CART_MAX];
187 FINT l_nx[CART_MAX], l_ny[CART_MAX], l_nz[CART_MAX];
188
189 CINTcart_comp(i_nx, i_ny, i_nz, i_l);
190 CINTcart_comp(j_nx, j_ny, j_nz, j_l);
191 CINTcart_comp(k_nx, k_ny, k_nz, k_l);
192 CINTcart_comp(l_nx, l_ny, l_nz, l_l);
193
194 ofx = 0;
195 ofy = envs->g_size;
196 ofz = envs->g_size * 2;
197 n = 0;
198 for (j = 0; j < nfj; j++) {
199 for (l = 0; l < nfl; l++) {
200 oflx = ofx + dj * j_nx[j] + dl * l_nx[l];
201 ofly = ofy + dj * j_ny[j] + dl * l_ny[l];
202 oflz = ofz + dj * j_nz[j] + dl * l_nz[l];
203 for (k = 0; k < nfk; k++) {
204 ofkx = oflx + dk * k_nx[k];
205 ofky = ofly + dk * k_ny[k];
206 ofkz = oflz + dk * k_nz[k];
207 switch (i_l) {
208 case 0:
209 idx[n+0] = ofkx;
210 idx[n+1] = ofky;
211 idx[n+2] = ofkz;
212 n += 3;
213 break;
214 case 1:
215 idx[n+0] = ofkx + di;
216 idx[n+1] = ofky;
217 idx[n+2] = ofkz;
218 idx[n+3] = ofkx;
219 idx[n+4] = ofky + di;
220 idx[n+5] = ofkz;
221 idx[n+6] = ofkx;
222 idx[n+7] = ofky;
223 idx[n+8] = ofkz + di;
224 n += 9;
225 break;
226 case 2:
227 idx[n+0 ] = ofkx + di*2;
228 idx[n+1 ] = ofky;
229 idx[n+2 ] = ofkz;
230 idx[n+3 ] = ofkx + di;
231 idx[n+4 ] = ofky + di;
232 idx[n+5 ] = ofkz;
233 idx[n+6 ] = ofkx + di;
234 idx[n+7 ] = ofky;
235 idx[n+8 ] = ofkz + di;
236 idx[n+9 ] = ofkx;
237 idx[n+10] = ofky + di*2;
238 idx[n+11] = ofkz;
239 idx[n+12] = ofkx;
240 idx[n+13] = ofky + di;
241 idx[n+14] = ofkz + di;
242 idx[n+15] = ofkx;
243 idx[n+16] = ofky;
244 idx[n+17] = ofkz + di*2;
245 n += 18;
246 break;
247 default:
248 for (i = 0; i < nfi; i++) {
249 idx[n+0] = ofkx + di * i_nx[i]; //(:,ix,kx,lx,jx,1)
250 idx[n+1] = ofky + di * i_ny[i]; //(:,iy,ky,ly,jy,2)
251 idx[n+2] = ofkz + di * i_nz[i]; //(:,iz,kz,lz,jz,3)
252 n += 3;
253 } // i
254 }
255 } // k
256 } // l
257 } // j
258 }
259
260
261 /*
262 * g(nroots,0:nmax,0:mmax)
263 */
CINTg0_2e_2d(double * g,struct _BC * bc,const CINTEnvVars * envs)264 void CINTg0_2e_2d(double *g, struct _BC *bc, const CINTEnvVars *envs)
265 {
266 const FINT nroots = envs->nrys_roots;
267 const FINT nmax = envs->li_ceil + envs->lj_ceil;
268 const FINT mmax = envs->lk_ceil + envs->ll_ceil;
269 const FINT dm = envs->g2d_klmax;
270 const FINT dn = envs->g2d_ijmax;
271 FINT i, j, m, n, off;
272 DEF_GXYZ(double, g, gx, gy, gz);
273 const double *c00;
274 const double *c0p;
275 const double *b01 = bc->b01;
276 const double *b00 = bc->b00;
277 const double *b10 = bc->b10;
278 double *p0x, *p0y, *p0z;
279 const double *p1x, *p1y, *p1z, *p2x, *p2y, *p2z;
280
281 for (i = 0; i < nroots; i++) {
282 gx[i] = 1;
283 gy[i] = 1;
284 //gz[i] = w[i];
285 }
286
287 if (nmax > 0) {
288 p0x = gx + dn;
289 p0y = gy + dn;
290 p0z = gz + dn;
291 p1x = gx - dn;
292 p1y = gy - dn;
293 p1z = gz - dn;
294 // gx(irys,0,1) = c00(irys) * gx(irys,0,0)
295 for (c00 = bc->c00, i = 0; i < nroots; i++, c00+=3) {
296 p0x[i] = c00[0] * gx[i];
297 p0y[i] = c00[1] * gy[i];
298 p0z[i] = c00[2] * gz[i];
299 }
300 // gx(irys,0,n+1) = c00(irys)*gx(irys,0,n)
301 // + n*b10(irys)*gx(irys,0,n-1)
302 for (n = 1; n < nmax; n++) {
303 off = n * dn;
304 for (c00 = bc->c00, i = 0, j = off;
305 i < nroots; i++, j++, c00+=3) {
306 p0x[j] = c00[0] * gx[j] + n * b10[i] * p1x[j];
307 p0y[j] = c00[1] * gy[j] + n * b10[i] * p1y[j];
308 p0z[j] = c00[2] * gz[j] + n * b10[i] * p1z[j];
309 }
310 }
311 }
312
313 if (mmax > 0) {
314 p0x = gx + dm;
315 p0y = gy + dm;
316 p0z = gz + dm;
317 p1x = gx - dm;
318 p1y = gy - dm;
319 p1z = gz - dm;
320 // gx(irys,1,0) = c0p(irys) * gx(irys,0,0)
321 for (c0p = bc->c0p, i = 0; i < nroots; i++, c0p+=3) {
322 p0x[i] = c0p[0] * gx[i];
323 p0y[i] = c0p[1] * gy[i];
324 p0z[i] = c0p[2] * gz[i];
325 }
326 // gx(irys,m+1,0) = c0p(irys)*gx(irys,m,0)
327 // + m*b01(irys)*gx(irys,m-1,0)
328 for (m = 1; m < mmax; m++) {
329 off = m * dm;
330 for (c0p = bc->c0p, i = 0, j = off;
331 i < nroots; i++, j++, c0p+=3) {
332 p0x[j] = c0p[0] * gx[j] + m * b01[i] * p1x[j];
333 p0y[j] = c0p[1] * gy[j] + m * b01[i] * p1y[j];
334 p0z[j] = c0p[2] * gz[j] + m * b01[i] * p1z[j];
335 }
336 }
337 }
338
339 if (nmax > 0 && mmax > 0) {
340 p0x = gx + dn;
341 p0y = gy + dn;
342 p0z = gz + dn;
343 p1x = gx - dn;
344 p1y = gy - dn;
345 p1z = gz - dn;
346 p2x = gx - dm;
347 p2y = gy - dm;
348 p2z = gz - dm;
349 // gx(irys,1,1) = c0p(irys)*gx(irys,0,1)
350 // + b00(irys)*gx(irys,0,0)
351 for (c0p = bc->c0p, i = 0; i < nroots; i++, c0p+=3) {
352 p0x[i+dm] = c0p[0] * p0x[i] + b00[i] * gx[i];
353 p0y[i+dm] = c0p[1] * p0y[i] + b00[i] * gy[i];
354 p0z[i+dm] = c0p[2] * p0z[i] + b00[i] * gz[i];
355 }
356
357 // gx(irys,m+1,1) = c0p(irys)*gx(irys,m,1)
358 // + m*b01(irys)*gx(irys,m-1,1)
359 // + b00(irys)*gx(irys,m,0)
360 for (m = 1; m < mmax; m++) {
361 off = m * dm + dn;
362 for (c0p = bc->c0p, i = 0, j = off;
363 i < nroots; i++, j++, c0p+=3) {
364 gx[j+dm] = c0p[0]*gx[j] + m*b01[i]*p2x[j] +b00[i]*p1x[j];
365 gy[j+dm] = c0p[1]*gy[j] + m*b01[i]*p2y[j] +b00[i]*p1y[j];
366 gz[j+dm] = c0p[2]*gz[j] + m*b01[i]*p2z[j] +b00[i]*p1z[j];
367 }
368 }
369
370 // gx(irys,m,n+1) = c00(irys)*gx(irys,m,n)
371 // + n*b10(irys)*gx(irys,m,n-1)
372 // + m*b00(irys)*gx(irys,m-1,n)
373 for (m = 1; m <= mmax; m++) {
374 for (n = 1; n < nmax; n++) {
375 off = m * dm + n * dn;
376 for (c00 = bc->c00, i = 0, j = off;
377 i < nroots; i++, j++, c00+=3) {
378 p0x[j] = c00[0]*gx[j] +n*b10[i]*p1x[j] + m*b00[i]*p2x[j];
379 p0y[j] = c00[1]*gy[j] +n*b10[i]*p1y[j] + m*b00[i]*p2y[j];
380 p0z[j] = c00[2]*gz[j] +n*b10[i]*p1z[j] + m*b00[i]*p2z[j];
381 }
382 }
383 }
384 }
385 }
386
387
388 /*
389 * g0[i,k,l,j] = < ik | lj > = ( i j | k l )
390 */
391 /* 2d is based on l,j */
CINTg0_lj2d_4d(double * g,const CINTEnvVars * envs)392 void CINTg0_lj2d_4d(double *g, const CINTEnvVars *envs)
393 {
394 const FINT nmax = envs->li_ceil + envs->lj_ceil;
395 const FINT mmax = envs->lk_ceil + envs->ll_ceil;
396 const FINT li = envs->li_ceil;
397 const FINT lk = envs->lk_ceil;
398 //const FINT ll = envs->ll_ceil;
399 const FINT lj = envs->lj_ceil;
400 const FINT nroots = envs->nrys_roots;
401 FINT i, j, k, l, ptr, n;
402 const FINT di = envs->g_stride_i;
403 const FINT dk = envs->g_stride_k;
404 const FINT dl = envs->g_stride_l;
405 const FINT dj = envs->g_stride_j;
406 const double *rirj = envs->rirj;
407 const double *rkrl = envs->rkrl;
408 DEF_GXYZ(double, g, gx, gy, gz);
409 const double *p1x, *p1y, *p1z, *p2x, *p2y, *p2z;
410 double rx, ry, rz;
411
412 // g(i,...,j) = rirj * g(i-1,...,j) + g(i-1,...,j+1)
413 rx = rirj[0];
414 ry = rirj[1];
415 rz = rirj[2];
416 p1x = gx - di;
417 p1y = gy - di;
418 p1z = gz - di;
419 p2x = gx - di + dj;
420 p2y = gy - di + dj;
421 p2z = gz - di + dj;
422 for (i = 1; i <= li; i++) {
423 for (j = 0; j <= nmax-i; j++) {
424 for (l = 0; l <= mmax; l++) {
425 ptr = j*dj + l*dl + i*di;
426 for (n = ptr; n < ptr+nroots; n++) {
427 gx[n] = rx * p1x[n] + p2x[n];
428 gy[n] = ry * p1y[n] + p2y[n];
429 gz[n] = rz * p1z[n] + p2z[n];
430 }
431 } } }
432
433 // g(...,k,l,..) = rkrl * g(...,k-1,l,..) + g(...,k-1,l+1,..)
434 rx = rkrl[0];
435 ry = rkrl[1];
436 rz = rkrl[2];
437 p1x = gx - dk;
438 p1y = gy - dk;
439 p1z = gz - dk;
440 p2x = gx - dk + dl;
441 p2y = gy - dk + dl;
442 p2z = gz - dk + dl;
443 for (j = 0; j <= lj; j++) {
444 for (k = 1; k <= lk; k++) {
445 for (l = 0; l <= mmax-k; l++) {
446 ptr = j*dj + l*dl + k*dk;
447 for (n = ptr; n < ptr+dk; n++) {
448 gx[n] = rx * p1x[n] + p2x[n];
449 gy[n] = ry * p1y[n] + p2y[n];
450 gz[n] = rz * p1z[n] + p2z[n];
451 }
452 } } }
453 }
454 /* 2d is based on k,j */
CINTg0_kj2d_4d(double * g,const CINTEnvVars * envs)455 void CINTg0_kj2d_4d(double *g, const CINTEnvVars *envs)
456 {
457 const FINT nmax = envs->li_ceil + envs->lj_ceil;
458 const FINT mmax = envs->lk_ceil + envs->ll_ceil;
459 const FINT li = envs->li_ceil;
460 //const FINT lk = envs->lk_ceil;
461 const FINT ll = envs->ll_ceil;
462 const FINT lj = envs->lj_ceil;
463 const FINT nroots = envs->nrys_roots;
464 FINT i, j, k, l, ptr, n;
465 const FINT di = envs->g_stride_i;
466 const FINT dk = envs->g_stride_k;
467 const FINT dl = envs->g_stride_l;
468 const FINT dj = envs->g_stride_j;
469 const double *rirj = envs->rirj;
470 const double *rkrl = envs->rkrl;
471 DEF_GXYZ(double, g, gx, gy, gz);
472 const double *p1x, *p1y, *p1z, *p2x, *p2y, *p2z;
473 double rx, ry, rz;
474
475 // g(i,...,j) = rirj * g(i-1,...,j) + g(i-1,...,j+1)
476 rx = rirj[0];
477 ry = rirj[1];
478 rz = rirj[2];
479 p1x = gx - di;
480 p1y = gy - di;
481 p1z = gz - di;
482 p2x = gx - di + dj;
483 p2y = gy - di + dj;
484 p2z = gz - di + dj;
485 for (i = 1; i <= li; i++) {
486 for (j = 0; j <= nmax-i; j++) {
487 for (k = 0; k <= mmax; k++) {
488 ptr = j*dj + k*dk + i*di;
489 for (n = ptr; n < ptr+nroots; n++) {
490 gx[n] = rx * p1x[n] + p2x[n];
491 gy[n] = ry * p1y[n] + p2y[n];
492 gz[n] = rz * p1z[n] + p2z[n];
493 }
494 } } }
495
496 // g(...,k,l,..) = rkrl * g(...,k,l-1,..) + g(...,k+1,l-1,..)
497 rx = rkrl[0];
498 ry = rkrl[1];
499 rz = rkrl[2];
500 p1x = gx - dl;
501 p1y = gy - dl;
502 p1z = gz - dl;
503 p2x = gx - dl + dk;
504 p2y = gy - dl + dk;
505 p2z = gz - dl + dk;
506 for (j = 0; j <= lj; j++) {
507 for (l = 1; l <= ll; l++) {
508 for (k = 0; k <= mmax-l; k++) {
509 ptr = j*dj + l*dl + k*dk;
510 for (n = ptr; n < ptr+dk; n++) {
511 gx[n] = rx * p1x[n] + p2x[n];
512 gy[n] = ry * p1y[n] + p2y[n];
513 gz[n] = rz * p1z[n] + p2z[n];
514 }
515 } } }
516 }
517 /* 2d is based on i,l */
CINTg0_il2d_4d(double * g,const CINTEnvVars * envs)518 void CINTg0_il2d_4d(double *g, const CINTEnvVars *envs)
519 {
520 const FINT nmax = envs->li_ceil + envs->lj_ceil;
521 const FINT mmax = envs->lk_ceil + envs->ll_ceil;
522 //const FINT li = envs->li_ceil;
523 const FINT lk = envs->lk_ceil;
524 const FINT ll = envs->ll_ceil;
525 const FINT lj = envs->lj_ceil;
526 const FINT nroots = envs->nrys_roots;
527 FINT i, j, k, l, ptr, n;
528 const FINT di = envs->g_stride_i;
529 const FINT dk = envs->g_stride_k;
530 const FINT dl = envs->g_stride_l;
531 const FINT dj = envs->g_stride_j;
532 const double *rirj = envs->rirj;
533 const double *rkrl = envs->rkrl;
534 DEF_GXYZ(double, g, gx, gy, gz);
535 const double *p1x, *p1y, *p1z, *p2x, *p2y, *p2z;
536 double rx, ry, rz;
537
538 // g(...,k,l,..) = rkrl * g(...,k-1,l,..) + g(...,k-1,l+1,..)
539 rx = rkrl[0];
540 ry = rkrl[1];
541 rz = rkrl[2];
542 p1x = gx - dk;
543 p1y = gy - dk;
544 p1z = gz - dk;
545 p2x = gx - dk + dl;
546 p2y = gy - dk + dl;
547 p2z = gz - dk + dl;
548 for (k = 1; k <= lk; k++) {
549 for (l = 0; l <= mmax-k; l++) {
550 for (i = 0; i <= nmax; i++) {
551 ptr = l*dl + k*dk + i*di;
552 for (n = ptr; n < ptr+nroots; n++) {
553 gx[n] = rx * p1x[n] + p2x[n];
554 gy[n] = ry * p1y[n] + p2y[n];
555 gz[n] = rz * p1z[n] + p2z[n];
556 }
557 } } }
558
559 // g(i,...,j) = rirj * g(i,...,j-1) + g(i+1,...,j-1)
560 rx = rirj[0];
561 ry = rirj[1];
562 rz = rirj[2];
563 p1x = gx - dj;
564 p1y = gy - dj;
565 p1z = gz - dj;
566 p2x = gx - dj + di;
567 p2y = gy - dj + di;
568 p2z = gz - dj + di;
569 for (j = 1; j <= lj; j++) {
570 for (l = 0; l <= ll; l++) {
571 for (k = 0; k <= lk; k++) {
572 ptr = j*dj + l*dl + k*dk;
573 for (n = ptr; n < ptr+dk-di*j; n++) {
574 gx[n] = rx * p1x[n] + p2x[n];
575 gy[n] = ry * p1y[n] + p2y[n];
576 gz[n] = rz * p1z[n] + p2z[n];
577 }
578 } } }
579 }
580 /* 2d is based on i,k */
CINTg0_ik2d_4d(double * g,const CINTEnvVars * envs)581 void CINTg0_ik2d_4d(double *g, const CINTEnvVars *envs)
582 {
583 const FINT nmax = envs->li_ceil + envs->lj_ceil;
584 const FINT mmax = envs->lk_ceil + envs->ll_ceil;
585 //const FINT li = envs->li_ceil;
586 const FINT lk = envs->lk_ceil;
587 const FINT ll = envs->ll_ceil;
588 const FINT lj = envs->lj_ceil;
589 const FINT nroots = envs->nrys_roots;
590 FINT i, j, k, l, ptr, n;
591 const FINT di = envs->g_stride_i;
592 const FINT dk = envs->g_stride_k;
593 const FINT dl = envs->g_stride_l;
594 const FINT dj = envs->g_stride_j;
595 const double *rirj = envs->rirj;
596 const double *rkrl = envs->rkrl;
597 DEF_GXYZ(double, g, gx, gy, gz);
598 const double *p1x, *p1y, *p1z, *p2x, *p2y, *p2z;
599 double rx, ry, rz;
600
601 // g(...,k,l,..) = rkrl * g(...,k,l-1,..) + g(...,k+1,l-1,..)
602 rx = rkrl[0];
603 ry = rkrl[1];
604 rz = rkrl[2];
605 p1x = gx - dl;
606 p1y = gy - dl;
607 p1z = gz - dl;
608 p2x = gx - dl + dk;
609 p2y = gy - dl + dk;
610 p2z = gz - dl + dk;
611 for (l = 1; l <= ll; l++) {
612 // (:,i) is full, so loop:k and loop:n can be merged to
613 // for(n = l*dl; n < ptr+dl-dk*l; n++)
614 for (k = 0; k <= mmax-l; k++) {
615 for (i = 0; i <= nmax; i++) {
616 ptr = l*dl + k*dk + i*di;
617 for (n = ptr; n < ptr+nroots; n++) {
618 gx[n] = rx * p1x[n] + p2x[n];
619 gy[n] = ry * p1y[n] + p2y[n];
620 gz[n] = rz * p1z[n] + p2z[n];
621 }
622 } }
623 }
624
625 // g(i,...,j) = rirj * g(i,...,j-1) + g(i+1,...,j-1)
626 rx = rirj[0];
627 ry = rirj[1];
628 rz = rirj[2];
629 p1x = gx - dj;
630 p1y = gy - dj;
631 p1z = gz - dj;
632 p2x = gx - dj + di;
633 p2y = gy - dj + di;
634 p2z = gz - dj + di;
635 for (j = 1; j <= lj; j++) {
636 for (l = 0; l <= ll; l++) {
637 for (k = 0; k <= lk; k++) {
638 ptr = j*dj + l*dl + k*dk;
639 for (n = ptr; n < ptr+dk-di*j; n++) {
640 gx[n] = rx * p1x[n] + p2x[n];
641 gy[n] = ry * p1y[n] + p2y[n];
642 gz[n] = rz * p1z[n] + p2z[n];
643 }
644 } } }
645 }
646 /************* some special g0_4d results *************/
647 /* 4 digits stand for i_ceil, k_ceil, l_ceil, j_ceil */
_g0_lj_4d_0001(double * g,double * c,const double * r)648 static inline void _g0_lj_4d_0001(double *g, double *c,
649 const double *r)
650 {
651 g[0] = 1;
652 g[1] = c[0];
653 g[2] = 1;
654 g[3] = c[1];
655 //g[4] = w[0];
656 g[5] = c[2] * g[4];
657 }
_g0_lj_4d_1000(double * g,double * c,const double * r)658 static inline void _g0_lj_4d_1000(double *g, double *c,
659 const double *r)
660 {
661 g[0] = 1;
662 g[1] = r[0] + c[0];
663 g[4] = 1;
664 g[5] = r[1] + c[1];
665 //g[8] = w[0];
666 g[9] =(r[2] + c[2]) * g[8];
667 }
_g0_lj_4d_0002(double * g,double * c,double * b,const double * r)668 static inline void _g0_lj_4d_0002(double *g, double *c, double *b,
669 const double *r)
670 {
671 g[0 ] = 1;
672 g[1 ] = 1;
673 g[2 ] = c[0];
674 g[3 ] = c[3];
675 g[4 ] = c[0] * c[0] + b[0];
676 g[5 ] = c[3] * c[3] + b[1];
677 g[6 ] = 1;
678 g[7 ] = 1;
679 g[8 ] = c[1];
680 g[9 ] = c[4];
681 g[10] = c[1] * c[1] + b[0];
682 g[11] = c[4] * c[4] + b[1];
683 //g[12] = w[0];
684 //g[13] = w[1];
685 g[14] = c[2] * g[12];
686 g[15] = c[5] * g[13];
687 g[16] =(c[2] * c[2] + b[0])* g[12];
688 g[17] =(c[5] * c[5] + b[1])* g[13];
689 }
_g0_lj_4d_1001(double * g,double * c,double * b,const double * r)690 static inline void _g0_lj_4d_1001(double *g, double *c, double *b,
691 const double *r)
692 {
693 double rc[] = {r[0]+c[0], r[0]+c[3],
694 r[1]+c[1], r[1]+c[4],
695 r[2]+c[2], r[2]+c[5]};
696 g[0 ] = 1;
697 g[1 ] = 1;
698 g[2 ] = rc[0];
699 g[3 ] = rc[1];
700 g[4 ] = c[0];
701 g[5 ] = c[3];
702 g[6 ] = rc[0] * c[0] + b[0];
703 g[7 ] = rc[1] * c[3] + b[1];
704 g[12] = 1;
705 g[13] = 1;
706 g[14] = rc[2];
707 g[15] = rc[3];
708 g[16] = c[1];
709 g[17] = c[4];
710 g[18] = rc[2] * c[1] + b[0];
711 g[19] = rc[3] * c[4] + b[1];
712 //g[24] = w[0];
713 //g[25] = w[1];
714 g[26] = rc[4] * g[24];
715 g[27] = rc[5] * g[25];
716 g[28] = c[2] * g[24];
717 g[29] = c[5] * g[25];
718 g[30] =(rc[4] * c[2] + b[0])* g[24];
719 g[31] =(rc[5] * c[5] + b[1])* g[25];
720 }
_g0_lj_4d_2000(double * g,double * c,double * b,const double * r)721 static inline void _g0_lj_4d_2000(double *g, double *c, double *b,
722 const double *r)
723 {
724 double rc[] = {r[0]+c[0], r[0]+c[3],
725 r[1]+c[1], r[1]+c[4],
726 r[2]+c[2], r[2]+c[5]};
727 g[0 ] = 1;
728 g[1 ] = 1;
729 g[2 ] = rc[0];
730 g[3 ] = rc[1];
731 g[4 ] = rc[0] * rc[0] + b[0];
732 g[5 ] = rc[1] * rc[1] + b[1];
733 g[18] = 1;
734 g[19] = 1;
735 g[20] = rc[2];
736 g[21] = rc[3];
737 g[22] = rc[2] * rc[2] + b[0];
738 g[23] = rc[3] * rc[3] + b[1];
739 //g[36] = w[0];
740 //g[37] = w[1];
741 g[38] = rc[4] * g[36];
742 g[39] = rc[5] * g[37];
743 g[40] =(rc[4] * rc[4] + b[0])* g[36];
744 g[41] =(rc[5] * rc[5] + b[1])* g[37];
745 }
_g0_lj_4d_0003(double * g,double * c,double * b,const double * r)746 static inline void _g0_lj_4d_0003(double *g, double *c, double *b,
747 const double *r)
748 {
749 g[0 ] = 1;
750 g[1 ] = 1;
751 g[2 ] = c[0];
752 g[3 ] = c[3];
753 g[4 ] = c[0] * c[0] + b[0];
754 g[5 ] = c[3] * c[3] + b[1];
755 g[6 ] = c[0] *(c[0] * c[0] + 3 * b[0]);
756 g[7 ] = c[3] *(c[3] * c[3] + 3 * b[1]);
757 g[8 ] = 1;
758 g[9 ] = 1;
759 g[10] = c[1];
760 g[11] = c[4];
761 g[12] = c[1] * c[1] + b[0];
762 g[13] = c[4] * c[4] + b[1];
763 g[14] = c[1] *(c[1] * c[1] + 3 * b[0]);
764 g[15] = c[4] *(c[4] * c[4] + 3 * b[1]);
765 //g[16] = w[0];
766 //g[17] = w[1];
767 g[18] = c[2] * g[16];
768 g[19] = c[5] * g[17];
769 g[20] =(c[2] * c[2] + b[0])* g[16];
770 g[21] =(c[5] * c[5] + b[1])* g[17];
771 g[22] =(c[2] * c[2] + 3 * b[0])* c[2] * g[16];
772 g[23] =(c[5] * c[5] + 3 * b[1])* c[5] * g[17];
773 }
_g0_lj_4d_1002(double * g,double * c,double * b,const double * r)774 static inline void _g0_lj_4d_1002(double *g, double *c, double *b,
775 const double *r)
776 {
777 double rc[] = {r[0]+c[0], r[0]+c[3],
778 r[1]+c[1], r[1]+c[4],
779 r[2]+c[2], r[2]+c[5]};
780 g[0 ] = 1;
781 g[1 ] = 1;
782 g[2 ] = rc[0];
783 g[3 ] = rc[1];
784 g[4 ] = c[0];
785 g[5 ] = c[3];
786 g[6 ] = rc[0] * c[0] + b[0];
787 g[7 ] = rc[1] * c[3] + b[1];
788 g[8 ] = c[0] * c[0] + b[0];
789 g[9 ] = c[3] * c[3] + b[1];
790 g[10] = rc[0]*c[0]*c[0] + b[0]*(rc[0]+2*c[0]);
791 g[11] = rc[1]*c[3]*c[3] + b[1]*(rc[1]+2*c[3]);
792 g[16] = 1;
793 g[17] = 1;
794 g[18] = rc[2];
795 g[19] = rc[3];
796 g[20] = c[1];
797 g[21] = c[4];
798 g[22] = rc[2] * c[1] + b[0];
799 g[23] = rc[3] * c[4] + b[1];
800 g[24] = c[1] * c[1] + b[0];
801 g[25] = c[4] * c[4] + b[1];
802 g[26] = rc[2]*c[1]*c[1] + b[0]*(rc[2]+2*c[1]);
803 g[27] = rc[3]*c[4]*c[4] + b[1]*(rc[3]+2*c[4]);
804 //g[32] = w[0];
805 //g[33] = w[1];
806 g[34] = rc[4] * g[32];
807 g[35] = rc[5] * g[33];
808 g[36] = c[2] * g[32];
809 g[37] = c[5] * g[33];
810 g[38] =(rc[4] * c[2] + b[0])* g[32];
811 g[39] =(rc[5] * c[5] + b[1])* g[33];
812 g[40] =(c[2] * c[2] + b[0])* g[32];
813 g[41] =(c[5] * c[5] + b[1])* g[33];
814 g[42] =(rc[4]*c[2]*c[2]+b[0]*(rc[4]+2*c[2]))*g[32];
815 g[43] =(rc[5]*c[5]*c[5]+b[1]*(rc[5]+2*c[5]))*g[33];
816 }
_g0_lj_4d_2001(double * g,double * c,double * b,const double * r)817 static inline void _g0_lj_4d_2001(double *g, double *c, double *b,
818 const double *r)
819 {
820 double rc[] = {r[0]+c[0], r[0]+c[3],
821 r[1]+c[1], r[1]+c[4],
822 r[2]+c[2], r[2]+c[5]};
823 g[0 ] = 1;
824 g[1 ] = 1;
825 g[2 ] = rc[0];
826 g[3 ] = rc[1];
827 g[4 ] = rc[0] * rc[0] + b[0];
828 g[5 ] = rc[1] * rc[1] + b[1];
829 g[6 ] = c[0];
830 g[7 ] = c[3];
831 g[8 ] = c[0] * rc[0] + b[0];
832 g[9 ] = c[3] * rc[1] + b[1];
833 g[10] = c[0]*rc[0]*rc[0] + b[0]*(2*rc[0]+c[0]);
834 g[11] = c[3]*rc[1]*rc[1] + b[1]*(2*rc[1]+c[3]);
835 g[24] = 1;
836 g[25] = 1;
837 g[26] = rc[2];
838 g[27] = rc[3];
839 g[28] = rc[2] * rc[2] + b[0];
840 g[29] = rc[3] * rc[3] + b[1];
841 g[30] = c[1];
842 g[31] = c[4];
843 g[32] = c[1] * rc[2] + b[0];
844 g[33] = c[4] * rc[3] + b[1];
845 g[34] = c[1]*rc[2]*rc[2] + b[0]*(2*rc[2]+c[1]);
846 g[35] = c[4]*rc[3]*rc[3] + b[1]*(2*rc[3]+c[4]);
847 //g[48] = w[0];
848 //g[49] = w[1];
849 g[50] = rc[4] * g[48];
850 g[51] = rc[5] * g[49];
851 g[52] =(rc[4] * rc[4] + b[0])* g[48];
852 g[53] =(rc[5] * rc[5] + b[1])* g[49];
853 g[54] = c[2] * g[48];
854 g[55] = c[5] * g[49];
855 g[56] =(c[2] * rc[4] + b[0])* g[48];
856 g[57] =(c[5] * rc[5] + b[1])* g[49];
857 g[58] =(c[2]*rc[4]*rc[4] + b[0]*(2*rc[4]+c[2]))* g[48];
858 g[59] =(c[5]*rc[5]*rc[5] + b[1]*(2*rc[5]+c[5]))* g[49];
859 }
_g0_lj_4d_3000(double * g,double * c,double * b,const double * r)860 static inline void _g0_lj_4d_3000(double *g, double *c, double *b,
861 const double *r)
862 {
863 double rc[] = {r[0]+c[0], r[0]+c[3],
864 r[1]+c[1], r[1]+c[4],
865 r[2]+c[2], r[2]+c[5]};
866 g[0 ] = 1;
867 g[1 ] = 1;
868 g[2 ] = rc[0];
869 g[3 ] = rc[1];
870 g[4 ] = rc[0] * rc[0] + b[0];
871 g[5 ] = rc[1] * rc[1] + b[1];
872 g[6 ] = rc[0] *(rc[0] * rc[0] + 3 * b[0]);
873 g[7 ] = rc[1] *(rc[1] * rc[1] + 3 * b[1]);
874 g[32] = 1;
875 g[33] = 1;
876 g[34] = rc[2];
877 g[35] = rc[3];
878 g[36] = rc[2] * rc[2] + b[0];
879 g[37] = rc[3] * rc[3] + b[1];
880 g[38] = rc[2] *(rc[2] * rc[2] + 3 * b[0]);
881 g[39] = rc[3] *(rc[3] * rc[3] + 3 * b[1]);
882 //g[64] = w[0];
883 //g[65] = w[1];
884 g[66] = rc[4] * g[64];
885 g[67] = rc[5] * g[65];
886 g[68] =(rc[4] * rc[4] + b[0])* g[64];
887 g[69] =(rc[5] * rc[5] + b[1])* g[65];
888 g[70] =(rc[4] * rc[4] + 3 * b[0])* rc[4] * g[64];
889 g[71] =(rc[5] * rc[5] + 3 * b[1])* rc[5] * g[65];
890 }
_g0_lj_4d_0011(double * g,double * c0,double * cp,double * b,const double * r0,const double * rp)891 static inline void _g0_lj_4d_0011(double *g, double *c0, double *cp, double *b,
892 const double *r0, const double *rp)
893 {
894 g[0 ] = 1;
895 g[1 ] = 1;
896 g[2 ] = cp[0];
897 g[3 ] = cp[3];
898 g[4 ] = c0[0];
899 g[5 ] = c0[3];
900 g[6 ] = cp[0] * c0[0] + b[0];
901 g[7 ] = cp[3] * c0[3] + b[1];
902 g[8 ] = 1;
903 g[9 ] = 1;
904 g[10] = cp[1];
905 g[11] = cp[4];
906 g[12] = c0[1];
907 g[13] = c0[4];
908 g[14] = cp[1] * c0[1] + b[0];
909 g[15] = cp[4] * c0[4] + b[1];
910 //g[16] = w[0];
911 //g[17] = w[1];
912 g[18] = cp[2] * g[16];
913 g[19] = cp[5] * g[17];
914 g[20] = c0[2] * g[16];
915 g[21] = c0[5] * g[17];
916 g[22] =(cp[2] * c0[2] + b[0]) * g[16];
917 g[23] =(cp[5] * c0[5] + b[1]) * g[17];
918 }
_g0_lj_4d_1010(double * g,double * c0,double * cp,double * b,const double * r0,const double * rp)919 static inline void _g0_lj_4d_1010(double *g, double *c0, double *cp, double *b,
920 const double *r0, const double *rp)
921 {
922 double rc[] = {r0[0]+c0[0], r0[0]+c0[3],
923 r0[1]+c0[1], r0[1]+c0[4],
924 r0[2]+c0[2], r0[2]+c0[5]};
925 g[0 ] = 1;
926 g[1 ] = 1;
927 g[2 ] = rc[0];
928 g[3 ] = rc[1];
929 g[4 ] = cp[0];
930 g[5 ] = cp[3];
931 g[6 ] = rc[0] * cp[0] + b[0];
932 g[7 ] = rc[1] * cp[3] + b[1];
933 g[16] = 1;
934 g[17] = 1;
935 g[18] = rc[2];
936 g[19] = rc[3];
937 g[20] = cp[1];
938 g[21] = cp[4];
939 g[22] = rc[2] * cp[1] + b[0];
940 g[23] = rc[3] * cp[4] + b[1];
941 //g[32] = w[0];
942 //g[33] = w[1];
943 g[34] = rc[4] * g[32];
944 g[35] = rc[5] * g[33];
945 g[36] = cp[2] * g[32];
946 g[37] = cp[5] * g[33];
947 g[38] =(rc[4]*cp[2] + b[0]) * g[32];
948 g[39] =(rc[5]*cp[5] + b[1]) * g[33];
949 }
_g0_lj_4d_0101(double * g,double * c0,double * cp,double * b,const double * r0,const double * rp)950 static inline void _g0_lj_4d_0101(double *g, double *c0, double *cp, double *b,
951 const double *r0, const double *rp)
952 {
953 double rc[] = {rp[0]+cp[0], rp[0]+cp[3],
954 rp[1]+cp[1], rp[1]+cp[4],
955 rp[2]+cp[2], rp[2]+cp[5]};
956 g[0 ] = 1;
957 g[1 ] = 1;
958 g[2 ] = rc[0];
959 g[3 ] = rc[1];
960 g[8 ] = c0[0];
961 g[9 ] = c0[3];
962 g[10] = rc[0] * c0[0] + b[0];
963 g[11] = rc[1] * c0[3] + b[1];
964 g[16] = 1;
965 g[17] = 1;
966 g[18] = rc[2];
967 g[19] = rc[3];
968 g[24] = c0[1];
969 g[25] = c0[4];
970 g[26] = rc[2] * c0[1] + b[0];
971 g[27] = rc[3] * c0[4] + b[1];
972 //g[32] = w[0];
973 //g[33] = w[1];
974 g[34] = rc[4] * g[32];
975 g[35] = rc[5] * g[33];
976 g[40] = c0[2] * g[32];
977 g[41] = c0[5] * g[33];
978 g[42] =(rc[4]*c0[2] + b[0]) * g[32];
979 g[43] =(rc[5]*c0[5] + b[1]) * g[33];
980 }
_g0_lj_4d_1100(double * g,double * c0,double * cp,double * b,const double * r0,const double * rp)981 static inline void _g0_lj_4d_1100(double *g, double *c0, double *cp, double *b,
982 const double *r0, const double *rp)
983 {
984 double rc0[] = {r0[0]+c0[0], r0[0]+c0[3],
985 r0[1]+c0[1], r0[1]+c0[4],
986 r0[2]+c0[2], r0[2]+c0[5]};
987 double rcp[] = {rp[0]+cp[0], rp[0]+cp[3],
988 rp[1]+cp[1], rp[1]+cp[4],
989 rp[2]+cp[2], rp[2]+cp[5]};
990 g[0 ] = 1;
991 g[1 ] = 1;
992 g[2 ] = rc0[0];
993 g[3 ] = rc0[1];
994 g[4 ] = rcp[0];
995 g[5 ] = rcp[1];
996 g[6 ] = rc0[0] * rcp[0] + b[0];
997 g[7 ] = rc0[1] * rcp[1] + b[1];
998 g[32] = 1;
999 g[33] = 1;
1000 g[34] = rc0[2];
1001 g[35] = rc0[3];
1002 g[36] = rcp[2];
1003 g[37] = rcp[3];
1004 g[38] = rc0[2] * rcp[2] + b[0];
1005 g[39] = rc0[3] * rcp[3] + b[1];
1006 //g[64] = w[0];
1007 //g[65] = w[1];
1008 g[66] = rc0[4] * g[64];
1009 g[67] = rc0[5] * g[65];
1010 g[68] = rcp[4] * g[64];
1011 g[69] = rcp[5] * g[65];
1012 g[70] =(rc0[4]*rcp[4] + b[0]) * g[64];
1013 g[71] =(rc0[5]*rcp[5] + b[1]) * g[65];
1014 }
_g0_lj_4d_0021(double * g,double * c0,double * cp,double * b0,double * b1)1015 static inline void _g0_lj_4d_0021(double *g, double *c0, double *cp,
1016 double *b0, double *b1)
1017 {
1018 g[0 ] = 1;
1019 g[1 ] = 1;
1020 g[2 ] = cp[0];
1021 g[3 ] = cp[3];
1022 g[4 ] = cp[0] * cp[0] + b1[0];
1023 g[5 ] = cp[3] * cp[3] + b1[1];
1024 g[6 ] = c0[0];
1025 g[7 ] = c0[3];
1026 g[8 ] = cp[0] * c0[0] + b0[0];
1027 g[9 ] = cp[3] * c0[3] + b0[1];
1028 g[10] = c0[0] * g[4] + 2 * b0[0] * cp[0];
1029 g[11] = c0[3] * g[5] + 2 * b0[1] * cp[3];
1030 g[12] = 1;
1031 g[13] = 1;
1032 g[14] = cp[1];
1033 g[15] = cp[4];
1034 g[16] = cp[1] * cp[1] + b1[0];
1035 g[17] = cp[4] * cp[4] + b1[1];
1036 g[18] = c0[1];
1037 g[19] = c0[4];
1038 g[20] = cp[1] * c0[1] + b0[0];
1039 g[21] = cp[4] * c0[4] + b0[1];
1040 g[22] = c0[1] * g[16] + 2 * b0[0] * cp[1];
1041 g[23] = c0[4] * g[17] + 2 * b0[1] * cp[4];
1042 //g[24] = w[0];
1043 //g[25] = w[1];
1044 g[26] = cp[2] * g[24];
1045 g[27] = cp[5] * g[25];
1046 g[28] =(cp[2] * cp[2] + b1[0]) * g[24];
1047 g[29] =(cp[5] * cp[5] + b1[1]) * g[25];
1048 g[30] = c0[2] * g[24];
1049 g[31] = c0[5] * g[25];
1050 g[32] =(cp[2] * c0[2] + b0[0]) * g[24];
1051 g[33] =(cp[5] * c0[5] + b0[1]) * g[25];
1052 g[34] = c0[2] * g[28] + 2 * b0[0] * g[26];
1053 g[35] = c0[5] * g[29] + 2 * b0[1] * g[27];
1054 }
_g0_lj_4d_1020(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1055 static inline void _g0_lj_4d_1020(double *g, double *c0, double *cp,
1056 double *b0, double *b1,
1057 const double *r0, const double *rp)
1058 {
1059 double rc[] = {r0[0]+c0[0], r0[0]+c0[3],
1060 r0[1]+c0[1], r0[1]+c0[4],
1061 r0[2]+c0[2], r0[2]+c0[5]};
1062 g[0 ] = 1;
1063 g[1 ] = 1;
1064 g[2 ] = rc[0];
1065 g[3 ] = rc[1];
1066 g[4 ] = cp[0];
1067 g[5 ] = cp[3];
1068 g[6 ] = cp[0] * rc[0] + b0[0];
1069 g[7 ] = cp[3] * rc[1] + b0[1];
1070 g[8 ] = cp[0] * cp[0] + b1[0];
1071 g[9 ] = cp[3] * cp[3] + b1[1];
1072 g[10] = rc[0] * g[8] + 2 * b0[0] * cp[0];
1073 g[11] = rc[1] * g[9] + 2 * b0[1] * cp[3];
1074 g[24] = 1;
1075 g[25] = 1;
1076 g[26] = rc[2];
1077 g[27] = rc[3];
1078 g[28] = cp[1];
1079 g[29] = cp[4];
1080 g[30] = cp[1] * rc[2] + b0[0];
1081 g[31] = cp[4] * rc[3] + b0[1];
1082 g[32] = cp[1] * cp[1] + b1[0];
1083 g[33] = cp[4] * cp[4] + b1[1];
1084 g[34] = rc[2] * g[32] + 2 * b0[0] * cp[1];
1085 g[35] = rc[3] * g[33] + 2 * b0[1] * cp[4];
1086 //g[48] = w[0];
1087 //g[49] = w[1];
1088 g[50] = rc[4] * g[48];
1089 g[51] = rc[5] * g[49];
1090 g[52] = cp[2] * g[48];
1091 g[53] = cp[5] * g[49];
1092 g[54] =(cp[2] * rc[4] + b0[0]) * g[48];
1093 g[55] =(cp[5] * rc[5] + b0[1]) * g[49];
1094 g[56] =(cp[2] * cp[2] + b1[0]) * g[48];
1095 g[57] =(cp[5] * cp[5] + b1[1]) * g[49];
1096 g[58] = rc[4] * g[56] + 2 * b0[0] * g[52];
1097 g[59] = rc[5] * g[57] + 2 * b0[1] * g[53];
1098 }
_g0_lj_4d_0111(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1099 static inline void _g0_lj_4d_0111(double *g, double *c0, double *cp,
1100 double *b0, double *b1,
1101 const double *r0, const double *rp)
1102 {
1103 double rc[] = {rp[0]+cp[0], rp[0]+cp[3],
1104 rp[1]+cp[1], rp[1]+cp[4],
1105 rp[2]+cp[2], rp[2]+cp[5]};
1106 g[0 ] = 1;
1107 g[1 ] = 1;
1108 g[2 ] = rc[0];
1109 g[3 ] = rc[1];
1110 g[4 ] = cp[0];
1111 g[5 ] = cp[3];
1112 g[6 ] = cp[0] * rc[0] + b1[0];
1113 g[7 ] = cp[3] * rc[1] + b1[1];
1114 g[12] = c0[0];
1115 g[13] = c0[3];
1116 g[14] = c0[0] * rc[0] + b0[0];
1117 g[15] = c0[3] * rc[1] + b0[1];
1118 g[16] = c0[0] * cp[0] + b0[0];
1119 g[17] = c0[3] * cp[3] + b0[1];
1120 g[18] = c0[0] * g[6] + b0[0] *(rc[0] + cp[0]);
1121 g[19] = c0[3] * g[7] + b0[1] *(rc[1] + cp[3]);
1122 g[24] = 1;
1123 g[25] = 1;
1124 g[26] = rc[2];
1125 g[27] = rc[3];
1126 g[28] = cp[1];
1127 g[29] = cp[4];
1128 g[30] = cp[1] * rc[2] + b1[0];
1129 g[31] = cp[4] * rc[3] + b1[1];
1130 g[36] = c0[1];
1131 g[37] = c0[4];
1132 g[38] = c0[1] * rc[2] + b0[0];
1133 g[39] = c0[4] * rc[3] + b0[1];
1134 g[40] = c0[1] * cp[1] + b0[0];
1135 g[41] = c0[4] * cp[4] + b0[1];
1136 g[42] = c0[1] * g[30] + b0[0] *(rc[2] + cp[1]);
1137 g[43] = c0[4] * g[31] + b0[1] *(rc[3] + cp[4]);
1138 //g[48] = w[0];
1139 //g[49] = w[1];
1140 g[50] = rc[4] * g[48];
1141 g[51] = rc[5] * g[49];
1142 g[52] = cp[2] * g[48];
1143 g[53] = cp[5] * g[49];
1144 g[54] =(cp[2] * rc[4] + b1[0]) * g[48];
1145 g[55] =(cp[5] * rc[5] + b1[1]) * g[49];
1146 g[60] = c0[2] * g[48];
1147 g[61] = c0[5] * g[49];
1148 g[62] =(c0[2] * rc[4] + b0[0]) * g[48];
1149 g[63] =(c0[5] * rc[5] + b0[1]) * g[49];
1150 g[64] =(c0[2] * cp[2] + b0[0]) * g[48];
1151 g[65] =(c0[5] * cp[5] + b0[1]) * g[49];
1152 g[66] = c0[2] * g[54] + b0[0] *(g[50] + g[52]);
1153 g[67] = c0[5] * g[55] + b0[1] *(g[51] + g[53]);
1154 }
_g0_lj_4d_1110(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1155 static inline void _g0_lj_4d_1110(double *g, double *c0, double *cp,
1156 double *b0, double *b1,
1157 const double *r0, const double *rp)
1158 {
1159 double rc0[] = {r0[0]+c0[0], r0[0]+c0[3],
1160 r0[1]+c0[1], r0[1]+c0[4],
1161 r0[2]+c0[2], r0[2]+c0[5]};
1162 double rcp[] = {rp[0]+cp[0], rp[0]+cp[3],
1163 rp[1]+cp[1], rp[1]+cp[4],
1164 rp[2]+cp[2], rp[2]+cp[5]};
1165 g[0 ] = 1;
1166 g[1 ] = 1;
1167 g[2 ] = rc0[0];
1168 g[3 ] = rc0[1];
1169 g[4 ] = rcp[0];
1170 g[5 ] = rcp[1];
1171 g[6 ] = rcp[0] * rc0[0] + b0[0];
1172 g[7 ] = rcp[1] * rc0[1] + b0[1];
1173 g[8 ] = cp[0];
1174 g[9 ] = cp[3];
1175 g[10] = cp[0] * rc0[0] + b0[0];
1176 g[11] = cp[3] * rc0[1] + b0[1];
1177 g[12] = cp[0] * rcp[0] + b1[0];
1178 g[13] = cp[3] * rcp[1] + b1[1];
1179 g[14] = rc0[0] * g[12] + b0[0] *(rcp[0] + cp[0]);
1180 g[15] = rc0[1] * g[13] + b0[1] *(rcp[1] + cp[3]);
1181 g[48] = 1;
1182 g[49] = 1;
1183 g[50] = rc0[2];
1184 g[51] = rc0[3];
1185 g[52] = rcp[2];
1186 g[53] = rcp[3];
1187 g[54] = rcp[2] * rc0[2] + b0[0];
1188 g[55] = rcp[3] * rc0[3] + b0[1];
1189 g[56] = cp[1];
1190 g[57] = cp[4];
1191 g[58] = cp[1] * rc0[2] + b0[0];
1192 g[59] = cp[4] * rc0[3] + b0[1];
1193 g[60] = cp[1] * rcp[2] + b1[0];
1194 g[61] = cp[4] * rcp[3] + b1[1];
1195 g[62] = rc0[2] * g[60] + b0[0] *(rcp[2] + cp[1]);
1196 g[63] = rc0[3] * g[61] + b0[1] *(rcp[3] + cp[4]);
1197 //g[96 ] = w[0];
1198 //g[97 ] = w[1];
1199 g[98 ] = rc0[4] * g[96 ];
1200 g[99 ] = rc0[5] * g[97 ];
1201 g[100] = rcp[4] * g[96 ];
1202 g[101] = rcp[5] * g[97 ];
1203 g[102] =(rcp[4] * rc0[4] + b0[0])* g[96 ];
1204 g[103] =(rcp[5] * rc0[5] + b0[1])* g[97 ];
1205 g[104] = cp[2] * g[96 ];
1206 g[105] = cp[5] * g[97 ];
1207 g[106] =(cp[2] * rc0[4] + b0[0])* g[96 ];
1208 g[107] =(cp[5] * rc0[5] + b0[1])* g[97 ];
1209 g[108] =(cp[2] * rcp[4] + b1[0])* g[96 ];
1210 g[109] =(cp[5] * rcp[5] + b1[1])* g[97 ];
1211 g[110] = rc0[4] * g[108] + b0[0] *(g[100] + g[104]);
1212 g[111] = rc0[5] * g[109] + b0[1] *(g[101] + g[105]);
1213 }
_g0_lj_4d_0201(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1214 static inline void _g0_lj_4d_0201(double *g, double *c0, double *cp,
1215 double *b0, double *b1,
1216 const double *r0, const double *rp)
1217 {
1218 double rc[] = {rp[0]+cp[0], rp[0]+cp[3],
1219 rp[1]+cp[1], rp[1]+cp[4],
1220 rp[2]+cp[2], rp[2]+cp[5]};
1221 g[0 ] = 1;
1222 g[1 ] = 1;
1223 g[2 ] = rc[0];
1224 g[3 ] = rc[1];
1225 g[4 ] = rc[0] * rc[0] + b1[0];
1226 g[5 ] = rc[1] * rc[1] + b1[1];
1227 g[18] = c0[0];
1228 g[19] = c0[3];
1229 g[20] = rc[0] * c0[0] + b0[0];
1230 g[21] = rc[1] * c0[3] + b0[1];
1231 g[22] = c0[0] * g[4] + 2 * b0[0] * rc[0];
1232 g[23] = c0[3] * g[5] + 2 * b0[1] * rc[1];
1233 g[36] = 1;
1234 g[37] = 1;
1235 g[38] = rc[2];
1236 g[39] = rc[3];
1237 g[40] = rc[2] * rc[2] + b1[0];
1238 g[41] = rc[3] * rc[3] + b1[1];
1239 g[54] = c0[1];
1240 g[55] = c0[4];
1241 g[56] = rc[2] * c0[1] + b0[0];
1242 g[57] = rc[3] * c0[4] + b0[1];
1243 g[58] = c0[1] * g[40] + 2 * b0[0] * rc[2];
1244 g[59] = c0[4] * g[41] + 2 * b0[1] * rc[3];
1245 //g[72] = w[0];
1246 //g[73] = w[1];
1247 g[74] = rc[4] * g[72];
1248 g[75] = rc[5] * g[73];
1249 g[76] =(rc[4] * rc[4] + b1[0])* g[72];
1250 g[77] =(rc[5] * rc[5] + b1[1])* g[73];
1251 g[90] = c0[2] * g[72];
1252 g[91] = c0[5] * g[73];
1253 g[92] =(rc[4] * c0[2] + b0[0])* g[72];
1254 g[93] =(rc[5] * c0[5] + b0[1])* g[73];
1255 g[94] = c0[2] * g[76] + 2 * b0[0] * g[74];
1256 g[95] = c0[5] * g[77] + 2 * b0[1] * g[75];
1257 }
_g0_lj_4d_1200(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1258 static inline void _g0_lj_4d_1200(double *g, double *c0, double *cp,
1259 double *b0, double *b1,
1260 const double *r0, const double *rp)
1261 {
1262 double rc0[] = {r0[0]+c0[0], r0[0]+c0[3],
1263 r0[1]+c0[1], r0[1]+c0[4],
1264 r0[2]+c0[2], r0[2]+c0[5]};
1265 double rcp[] = {rp[0]+cp[0], rp[0]+cp[3],
1266 rp[1]+cp[1], rp[1]+cp[4],
1267 rp[2]+cp[2], rp[2]+cp[5]};
1268 g[0 ] = 1;
1269 g[1 ] = 1;
1270 g[2 ] = rc0[0];
1271 g[3 ] = rc0[1];
1272 g[4 ] = rcp[0];
1273 g[5 ] = rcp[1];
1274 g[6 ] = rcp[0] * rc0[0] + b0[0];
1275 g[7 ] = rcp[1] * rc0[1] + b0[1];
1276 g[8 ] = rcp[0] * rcp[0] + b1[0];
1277 g[9 ] = rcp[1] * rcp[1] + b1[1];
1278 g[10] = rc0[0] * g[8] + 2 * b0[0] * rcp[0];
1279 g[11] = rc0[1] * g[9] + 2 * b0[1] * rcp[1];
1280 g[72] = 1;
1281 g[73] = 1;
1282 g[74] = rc0[2];
1283 g[75] = rc0[3];
1284 g[76] = rcp[2];
1285 g[77] = rcp[3];
1286 g[78] = rcp[2] * rc0[2] + b0[0];
1287 g[79] = rcp[3] * rc0[3] + b0[1];
1288 g[80] = rcp[2] * rcp[2] + b1[0];
1289 g[81] = rcp[3] * rcp[3] + b1[1];
1290 g[82] = rc0[2] * g[80] + 2 * b0[0] * rcp[2];
1291 g[83] = rc0[3] * g[81] + 2 * b0[1] * rcp[3];
1292 //g[144] = w[0];
1293 //g[145] = w[1];
1294 g[146] = rc0[4] * g[144];
1295 g[147] = rc0[5] * g[145];
1296 g[148] = rcp[4] * g[144];
1297 g[149] = rcp[5] * g[145];
1298 g[150] =(rcp[4] * rc0[4] + b0[0])* g[144];
1299 g[151] =(rcp[5] * rc0[5] + b0[1])* g[145];
1300 g[152] =(rcp[4] * rcp[4] + b1[0])* g[144];
1301 g[153] =(rcp[5] * rcp[5] + b1[1])* g[145];
1302 g[154] = rc0[4] * g[152] + 2 * b0[0] * g[148];
1303 g[155] = rc0[5] * g[153] + 2 * b0[1] * g[149];
1304 }
_g0_lj_4d_0012(double * g,double * c0,double * cp,double * b0,double * b1)1305 static inline void _g0_lj_4d_0012(double *g, double *c0, double *cp,
1306 double *b0, double *b1)
1307 {
1308 g[0 ] = 1;
1309 g[1 ] = 1;
1310 g[2 ] = cp[0];
1311 g[3 ] = cp[3];
1312 g[4 ] = c0[0];
1313 g[5 ] = c0[3];
1314 g[6 ] = cp[0] * c0[0] + b0[0];
1315 g[7 ] = cp[3] * c0[3] + b0[1];
1316 g[8 ] = c0[0] * c0[0] + b1[0];
1317 g[9 ] = c0[3] * c0[3] + b1[1];
1318 g[10] = cp[0] * g[8] + 2 * b0[0] * c0[0];
1319 g[11] = cp[3] * g[9] + 2 * b0[1] * c0[3];
1320 g[12] = 1;
1321 g[13] = 1;
1322 g[14] = cp[1];
1323 g[15] = cp[4];
1324 g[16] = c0[1];
1325 g[17] = c0[4];
1326 g[18] = cp[1] * c0[1] + b0[0];
1327 g[19] = cp[4] * c0[4] + b0[1];
1328 g[20] = c0[1] * c0[1] + b1[0];
1329 g[21] = c0[4] * c0[4] + b1[1];
1330 g[22] = cp[1] * g[20] + 2 * b0[0] * c0[1];
1331 g[23] = cp[4] * g[21] + 2 * b0[1] * c0[4];
1332 //g[24] = w[0];
1333 //g[25] = w[1];
1334 g[26] = cp[2] * g[24];
1335 g[27] = cp[5] * g[25];
1336 g[28] = c0[2] * g[24];
1337 g[29] = c0[5] * g[25];
1338 g[30] =(cp[2] * c0[2] + b0[0]) * g[24];
1339 g[31] =(cp[5] * c0[5] + b0[1]) * g[25];
1340 g[32] =(c0[2] * c0[2] + b1[0]) * g[24];
1341 g[33] =(c0[5] * c0[5] + b1[1]) * g[25];
1342 g[34] = cp[2] * g[32] + 2 * b0[0] * g[28];
1343 g[35] = cp[5] * g[33] + 2 * b0[1] * g[29];
1344 }
_g0_lj_4d_1011(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1345 static inline void _g0_lj_4d_1011(double *g, double *c0, double *cp,
1346 double *b0, double *b1,
1347 const double *r0, const double *rp)
1348 {
1349 double rc[] = {r0[0]+c0[0], r0[0]+c0[3],
1350 r0[1]+c0[1], r0[1]+c0[4],
1351 r0[2]+c0[2], r0[2]+c0[5]};
1352 g[0 ] = 1;
1353 g[1 ] = 1;
1354 g[2 ] = rc[0];
1355 g[3 ] = rc[1];
1356 g[4 ] = cp[0];
1357 g[5 ] = cp[3];
1358 g[6 ] = cp[0] * rc[0] + b0[0];
1359 g[7 ] = cp[3] * rc[1] + b0[1];
1360 g[8 ] = c0[0];
1361 g[9 ] = c0[3];
1362 g[10] = c0[0] * rc[0] + b1[0];
1363 g[11] = c0[3] * rc[1] + b1[1];
1364 g[12] = c0[0] * cp[0] + b0[0];
1365 g[13] = c0[3] * cp[3] + b0[1];
1366 g[14] = cp[0] * g[10] + b0[0] *(rc[0] + c0[0]);
1367 g[15] = cp[3] * g[11] + b0[1] *(rc[1] + c0[3]);
1368 g[24] = 1;
1369 g[25] = 1;
1370 g[26] = rc[2];
1371 g[27] = rc[3];
1372 g[28] = cp[1];
1373 g[29] = cp[4];
1374 g[30] = cp[1] * rc[2] + b0[0];
1375 g[31] = cp[4] * rc[3] + b0[1];
1376 g[32] = c0[1];
1377 g[33] = c0[4];
1378 g[34] = c0[1] * rc[2] + b1[0];
1379 g[35] = c0[4] * rc[3] + b1[1];
1380 g[36] = c0[1] * cp[1] + b0[0];
1381 g[37] = c0[4] * cp[4] + b0[1];
1382 g[38] = cp[1] * g[34] + b0[0] *(rc[2] + c0[1]);
1383 g[39] = cp[4] * g[35] + b0[1] *(rc[3] + c0[4]);
1384 //g[48] = w[0];
1385 //g[49] = w[1];
1386 g[50] = rc[4] * g[48];
1387 g[51] = rc[5] * g[49];
1388 g[52] = cp[2] * g[48];
1389 g[53] = cp[5] * g[49];
1390 g[54] =(cp[2] * rc[4] + b0[0])* g[48];
1391 g[55] =(cp[5] * rc[5] + b0[1])* g[49];
1392 g[56] = c0[2] * g[48];
1393 g[57] = c0[5] * g[49];
1394 g[58] =(c0[2] * rc[4] + b1[0])* g[48];
1395 g[59] =(c0[5] * rc[5] + b1[1])* g[49];
1396 g[60] =(c0[2] * cp[2] + b0[0])* g[48];
1397 g[61] =(c0[5] * cp[5] + b0[1])* g[49];
1398 g[62] = cp[2] * g[58] + b0[0] *(g[50] + g[56]);
1399 g[63] = cp[5] * g[59] + b0[1] *(g[51] + g[57]);
1400 }
_g0_lj_4d_2010(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1401 static inline void _g0_lj_4d_2010(double *g, double *c0, double *cp,
1402 double *b0, double *b1,
1403 const double *r0, const double *rp)
1404 {
1405 double rc[] = {r0[0]+c0[0], r0[0]+c0[3],
1406 r0[1]+c0[1], r0[1]+c0[4],
1407 r0[2]+c0[2], r0[2]+c0[5]};
1408 g[0 ] = 1;
1409 g[1 ] = 1;
1410 g[2 ] = rc[0];
1411 g[3 ] = rc[1];
1412 g[4 ] = rc[0] * rc[0] + b1[0];
1413 g[5 ] = rc[1] * rc[1] + b1[1];
1414 g[6 ] = cp[0];
1415 g[7 ] = cp[3];
1416 g[8 ] = cp[0] * rc[0] + b0[0];
1417 g[9 ] = cp[3] * rc[1] + b0[1];
1418 g[10] = cp[0] * g[4] + 2 * b0[0] * rc[0];
1419 g[11] = cp[3] * g[5] + 2 * b0[1] * rc[1];
1420 g[36] = 1;
1421 g[37] = 1;
1422 g[38] = rc[2];
1423 g[39] = rc[3];
1424 g[40] = rc[2] * rc[2] + b1[0];
1425 g[41] = rc[3] * rc[3] + b1[1];
1426 g[42] = cp[1];
1427 g[43] = cp[4];
1428 g[44] = cp[1] * rc[2] + b0[0];
1429 g[45] = cp[4] * rc[3] + b0[1];
1430 g[46] = cp[1] * g[40] + 2 * b0[0] * rc[2];
1431 g[47] = cp[4] * g[41] + 2 * b0[1] * rc[3];
1432 //g[72] = w[0];
1433 //g[73] = w[1];
1434 g[74] = rc[4] * g[72];
1435 g[75] = rc[5] * g[73];
1436 g[76] =(rc[4] * rc[4] + b1[0]) * g[72];
1437 g[77] =(rc[5] * rc[5] + b1[1]) * g[73];
1438 g[78] = cp[2] * g[72];
1439 g[79] = cp[5] * g[73];
1440 g[80] =(cp[2] * rc[4] + b0[0]) * g[72];
1441 g[81] =(cp[5] * rc[5] + b0[1]) * g[73];
1442 g[82] = cp[2] * g[76] + 2 * b0[0] * g[74];
1443 g[83] = cp[5] * g[77] + 2 * b0[1] * g[75];
1444 }
_g0_lj_4d_0102(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1445 static inline void _g0_lj_4d_0102(double *g, double *c0, double *cp,
1446 double *b0, double *b1,
1447 const double *r0, const double *rp)
1448 {
1449 double rc[] = {rp[0]+cp[0], rp[0]+cp[3],
1450 rp[1]+cp[1], rp[1]+cp[4],
1451 rp[2]+cp[2], rp[2]+cp[5]};
1452 g[0 ] = 1;
1453 g[1 ] = 1;
1454 g[2 ] = rc[0];
1455 g[3 ] = rc[1];
1456 g[8 ] = c0[0];
1457 g[9 ] = c0[3];
1458 g[10] = rc[0] * c0[0] + b0[0];
1459 g[11] = rc[1] * c0[3] + b0[1];
1460 g[16] = c0[0] * c0[0] + b1[0];
1461 g[17] = c0[3] * c0[3] + b1[1];
1462 g[18] = rc[0] * g[16] + 2 * b0[0] * c0[0];
1463 g[19] = rc[1] * g[17] + 2 * b0[1] * c0[3];
1464 g[24] = 1;
1465 g[25] = 1;
1466 g[26] = rc[2];
1467 g[27] = rc[3];
1468 g[32] = c0[1];
1469 g[33] = c0[4];
1470 g[34] = rc[2] * c0[1] + b0[0];
1471 g[35] = rc[3] * c0[4] + b0[1];
1472 g[40] = c0[1] * c0[1] + b1[0];
1473 g[41] = c0[4] * c0[4] + b1[1];
1474 g[42] = rc[2] * g[40] + 2 * b0[0] * c0[1];
1475 g[43] = rc[3] * g[41] + 2 * b0[1] * c0[4];
1476 //g[48] = w[0];
1477 //g[49] = w[1];
1478 g[50] = rc[4] * g[48];
1479 g[51] = rc[5] * g[49];
1480 g[56] = c0[2] * g[48];
1481 g[57] = c0[5] * g[49];
1482 g[58] =(rc[4] * c0[2] + b0[0]) * g[48];
1483 g[59] =(rc[5] * c0[5] + b0[1]) * g[49];
1484 g[64] =(c0[2] * c0[2] + b1[0]) * g[48];
1485 g[65] =(c0[5] * c0[5] + b1[1]) * g[49];
1486 g[66] = rc[4] * g[64] + 2 * b0[0] * g[56];
1487 g[67] = rc[5] * g[65] + 2 * b0[1] * g[57];
1488 }
_g0_lj_4d_1101(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1489 static inline void _g0_lj_4d_1101(double *g, double *c0, double *cp,
1490 double *b0, double *b1,
1491 const double *r0, const double *rp)
1492 {
1493 double rc0[] = {r0[0]+c0[0], r0[0]+c0[3],
1494 r0[1]+c0[1], r0[1]+c0[4],
1495 r0[2]+c0[2], r0[2]+c0[5]};
1496 double rcp[] = {rp[0]+cp[0], rp[0]+cp[3],
1497 rp[1]+cp[1], rp[1]+cp[4],
1498 rp[2]+cp[2], rp[2]+cp[5]};
1499 g[0 ] = 1;
1500 g[1 ] = 1;
1501 g[2 ] = rc0[0];
1502 g[3 ] = rc0[1];
1503 g[4 ] = rcp[0];
1504 g[5 ] = rcp[1];
1505 g[6 ] = rcp[0] * rc0[0] + b0[0];
1506 g[7 ] = rcp[1] * rc0[1] + b0[1];
1507 g[16] = c0[0];
1508 g[17] = c0[3];
1509 g[18] = c0[0] * rc0[0] + b1[0];
1510 g[19] = c0[3] * rc0[1] + b1[1];
1511 g[20] = c0[0] * rcp[0] + b0[0];
1512 g[21] = c0[3] * rcp[1] + b0[1];
1513 g[22] = rcp[0] * g[18] + b0[0] *(rc0[0] + c0[0]);
1514 g[23] = rcp[1] * g[19] + b0[1] *(rc0[1] + c0[3]);
1515 g[48] = 1;
1516 g[49] = 1;
1517 g[50] = rc0[2];
1518 g[51] = rc0[3];
1519 g[52] = rcp[2];
1520 g[53] = rcp[3];
1521 g[54] = rcp[2] * rc0[2] + b0[0];
1522 g[55] = rcp[3] * rc0[3] + b0[1];
1523 g[64] = c0[1];
1524 g[65] = c0[4];
1525 g[66] = c0[1] * rc0[2] + b1[0];
1526 g[67] = c0[4] * rc0[3] + b1[1];
1527 g[68] = c0[1] * rcp[2] + b0[0];
1528 g[69] = c0[4] * rcp[3] + b0[1];
1529 g[70] = rcp[2] * g[66] + b0[0] *(rc0[2] + c0[1]);
1530 g[71] = rcp[3] * g[67] + b0[1] *(rc0[3] + c0[4]);
1531 //g[96 ] = w[0];
1532 //g[97 ] = w[1];
1533 g[98 ] = rc0[4] * g[96];
1534 g[99 ] = rc0[5] * g[97];
1535 g[100] = rcp[4] * g[96];
1536 g[101] = rcp[5] * g[97];
1537 g[102] =(rcp[4] * rc0[4] + b0[0])* g[96];
1538 g[103] =(rcp[5] * rc0[5] + b0[1])* g[97];
1539 g[112] = c0[2] * g[96];
1540 g[113] = c0[5] * g[97];
1541 g[114] =(c0[2] * rc0[4] + b1[0])* g[96];
1542 g[115] =(c0[5] * rc0[5] + b1[1])* g[97];
1543 g[116] =(c0[2] * rcp[4] + b0[0])* g[96];
1544 g[117] =(c0[5] * rcp[5] + b0[1])* g[97];
1545 g[118] = rcp[4] * g[114] + b0[0] *(g[98] + g[112]);
1546 g[119] = rcp[5] * g[115] + b0[1] *(g[99] + g[113]);
1547 }
_g0_lj_4d_2100(double * g,double * c0,double * cp,double * b0,double * b1,const double * r0,const double * rp)1548 static inline void _g0_lj_4d_2100(double *g, double *c0, double *cp,
1549 double *b0, double *b1,
1550 const double *r0, const double *rp)
1551 {
1552 double rc0[] = {r0[0]+c0[0], r0[0]+c0[3],
1553 r0[1]+c0[1], r0[1]+c0[4],
1554 r0[2]+c0[2], r0[2]+c0[5]};
1555 double rcp[] = {rp[0]+cp[0], rp[0]+cp[3],
1556 rp[1]+cp[1], rp[1]+cp[4],
1557 rp[2]+cp[2], rp[2]+cp[5]};
1558 g[0 ] = 1;
1559 g[1 ] = 1;
1560 g[2 ] = rc0[0];
1561 g[3 ] = rc0[1];
1562 g[4 ] = rc0[0] * rc0[0] + b1[0];
1563 g[5 ] = rc0[1] * rc0[1] + b1[1];
1564 g[6 ] = rcp[0];
1565 g[7 ] = rcp[1];
1566 g[8 ] = rcp[0] * rc0[0] + b0[0];
1567 g[9 ] = rcp[1] * rc0[1] + b0[1];
1568 g[10] = rcp[0] * g[4] + 2 * b0[0] * rc0[0];
1569 g[11] = rcp[1] * g[5] + 2 * b0[1] * rc0[1];
1570 g[72] = 1;
1571 g[73] = 1;
1572 g[74] = rc0[2];
1573 g[75] = rc0[3];
1574 g[76] = rc0[2] * rc0[2] + b1[0];
1575 g[77] = rc0[3] * rc0[3] + b1[1];
1576 g[78] = rcp[2];
1577 g[79] = rcp[3];
1578 g[80] = rcp[2] * rc0[2] + b0[0];
1579 g[81] = rcp[3] * rc0[3] + b0[1];
1580 g[82] = rcp[2] * g[76] + 2 * b0[0] * rc0[2];
1581 g[83] = rcp[3] * g[77] + 2 * b0[1] * rc0[3];
1582 //g[144] = w[0];
1583 //g[145] = w[1];
1584 g[146] = rc0[4] * g[144];
1585 g[147] = rc0[5] * g[145];
1586 g[148] =(rc0[4] * rc0[4] + b1[0])* g[144];
1587 g[149] =(rc0[5] * rc0[5] + b1[1])* g[145];
1588 g[150] = rcp[4] * g[144];
1589 g[151] = rcp[5] * g[145];
1590 g[152] =(rcp[4] * rc0[4] + b0[0])* g[144];
1591 g[153] =(rcp[5] * rc0[5] + b0[1])* g[145];
1592 g[154] = rcp[4] * g[148] + 2 * b0[0] * g[146];
1593 g[155] = rcp[5] * g[149] + 2 * b0[1] * g[147];
1594 }
1595 /************** end special g0_4d results *************/
1596
1597
1598
CINTg0_2e_lj2d4d(double * g,struct _BC * bc,const CINTEnvVars * envs)1599 void CINTg0_2e_lj2d4d(double *g, struct _BC *bc, const CINTEnvVars *envs)
1600 {
1601 const FINT nmax = envs->li_ceil + envs->lj_ceil;
1602 const FINT mmax = envs->lk_ceil + envs->ll_ceil;
1603 switch (nmax) {
1604 case 0: switch(mmax) {
1605 case 0: goto _g0_4d_default; // ssss
1606 case 1: switch (envs->lk_ceil) {
1607 case 0: _g0_lj_4d_0001(g, bc->c0p, envs->rkrl); goto normal_end;
1608 case 1: _g0_lj_4d_1000(g, bc->c0p, envs->rkrl); goto normal_end;
1609 default: goto error; }
1610 case 2: switch (envs->lk_ceil) {
1611 case 0: _g0_lj_4d_0002(g, bc->c0p, bc->b01, envs->rkrl); goto normal_end;
1612 case 1: _g0_lj_4d_1001(g, bc->c0p, bc->b01, envs->rkrl); goto normal_end;
1613 case 2: _g0_lj_4d_2000(g, bc->c0p, bc->b01, envs->rkrl); goto normal_end;
1614 default: goto error; }
1615 case 3: switch (envs->lk_ceil) {
1616 case 0: _g0_lj_4d_0003(g, bc->c0p, bc->b01, envs->rkrl); goto normal_end;
1617 case 1: _g0_lj_4d_1002(g, bc->c0p, bc->b01, envs->rkrl); goto normal_end;
1618 case 2: _g0_lj_4d_2001(g, bc->c0p, bc->b01, envs->rkrl); goto normal_end;
1619 case 3: _g0_lj_4d_3000(g, bc->c0p, bc->b01, envs->rkrl); goto normal_end;
1620 default: goto error; }
1621 default: goto _g0_4d_default; }
1622 case 1: switch(mmax) {
1623 case 0: switch (envs->li_ceil) {
1624 case 0: _g0_lj_4d_0001(g, bc->c00, envs->rirj); goto normal_end;
1625 case 1: _g0_lj_4d_1000(g, bc->c00, envs->rirj); goto normal_end;
1626 default: goto error; }
1627 case 1: switch (envs->lk_ceil) {
1628 case 0: switch (envs->li_ceil) {
1629 case 0: _g0_lj_4d_0011(g, bc->c00, bc->c0p, bc->b00, envs->rirj, envs->rkrl); goto normal_end;
1630 case 1: _g0_lj_4d_1010(g, bc->c00, bc->c0p, bc->b00, envs->rirj, envs->rkrl); goto normal_end;
1631 default: goto error; }
1632 case 1: switch (envs->li_ceil) {
1633 case 0: _g0_lj_4d_0101(g, bc->c00, bc->c0p, bc->b00, envs->rirj, envs->rkrl); goto normal_end;
1634 case 1: _g0_lj_4d_1100(g, bc->c00, bc->c0p, bc->b00, envs->rirj, envs->rkrl); goto normal_end;
1635 default: goto error; }
1636 default: goto error; }
1637 case 2: switch (envs->lk_ceil) {
1638 case 0: switch (envs->li_ceil) {
1639 case 0: _g0_lj_4d_0021(g, bc->c00, bc->c0p, bc->b00, bc->b01); goto normal_end;
1640 case 1: _g0_lj_4d_1020(g, bc->c00, bc->c0p, bc->b00, bc->b01, envs->rirj, envs->rkrl); goto normal_end;
1641 default: goto error; }
1642 case 1: switch (envs->li_ceil) {
1643 case 0: _g0_lj_4d_0111(g, bc->c00, bc->c0p, bc->b00, bc->b01, envs->rirj, envs->rkrl); goto normal_end;
1644 case 1: _g0_lj_4d_1110(g, bc->c00, bc->c0p, bc->b00, bc->b01, envs->rirj, envs->rkrl); goto normal_end;
1645 default: goto error; }
1646 case 2: switch (envs->li_ceil) {
1647 case 0: _g0_lj_4d_0201(g, bc->c00, bc->c0p, bc->b00, bc->b01, envs->rirj, envs->rkrl); goto normal_end;
1648 case 1: _g0_lj_4d_1200(g, bc->c00, bc->c0p, bc->b00, bc->b01, envs->rirj, envs->rkrl); goto normal_end;
1649 default: goto error; }
1650 default: goto error; }
1651 default: goto _g0_4d_default; }
1652 case 2: switch(mmax) {
1653 case 0: switch (envs->li_ceil) {
1654 case 0: _g0_lj_4d_0002(g, bc->c00, bc->b10, envs->rirj); goto normal_end;
1655 case 1: _g0_lj_4d_1001(g, bc->c00, bc->b10, envs->rirj); goto normal_end;
1656 case 2: _g0_lj_4d_2000(g, bc->c00, bc->b10, envs->rirj); goto normal_end;
1657 default: goto error; }
1658 case 1: switch (envs->lk_ceil) {
1659 case 0: switch (envs->li_ceil) {
1660 case 0: _g0_lj_4d_0012(g, bc->c00, bc->c0p, bc->b00, bc->b10); goto normal_end;
1661 case 1: _g0_lj_4d_1011(g, bc->c00, bc->c0p, bc->b00, bc->b10, envs->rirj, envs->rkrl); goto normal_end;
1662 case 2: _g0_lj_4d_2010(g, bc->c00, bc->c0p, bc->b00, bc->b10, envs->rirj, envs->rkrl); goto normal_end;
1663 default: goto error; }
1664 case 1: switch (envs->li_ceil) {
1665 case 0: _g0_lj_4d_0102(g, bc->c00, bc->c0p, bc->b00, bc->b10, envs->rirj, envs->rkrl); goto normal_end;
1666 case 1: _g0_lj_4d_1101(g, bc->c00, bc->c0p, bc->b00, bc->b10, envs->rirj, envs->rkrl); goto normal_end;
1667 case 2: _g0_lj_4d_2100(g, bc->c00, bc->c0p, bc->b00, bc->b10, envs->rirj, envs->rkrl); goto normal_end;
1668 default: goto error; }
1669 default: goto error; }
1670 default: goto _g0_4d_default; }
1671 case 3: switch(mmax) {
1672 case 0: switch (envs->li_ceil) {
1673 case 0: _g0_lj_4d_0003(g, bc->c00, bc->b10, envs->rirj); goto normal_end;
1674 case 1: _g0_lj_4d_1002(g, bc->c00, bc->b10, envs->rirj); goto normal_end;
1675 case 2: _g0_lj_4d_2001(g, bc->c00, bc->b10, envs->rirj); goto normal_end;
1676 case 3: _g0_lj_4d_3000(g, bc->c00, bc->b10, envs->rirj); goto normal_end;
1677 default: goto error; }
1678 default: goto _g0_4d_default; }
1679 default:
1680 _g0_4d_default:
1681 CINTg0_2e_2d(g, bc, envs);
1682 CINTg0_lj2d_4d(g, envs);
1683 }
1684 normal_end:
1685 return;
1686 error:
1687 fprintf(stderr, "Dimension error for CINTg0_2e_lj2d4d: iklj = %d %d %d %d",
1688 (int)envs->li_ceil, (int)envs->lk_ceil,
1689 (int)envs->ll_ceil, (int)envs->lj_ceil);
1690 #ifdef DEBUG
1691 exit(1);
1692 #endif
1693 }
1694
CINTg0_2e_kj2d4d(double * g,struct _BC * bc,const CINTEnvVars * envs)1695 void CINTg0_2e_kj2d4d(double *g, struct _BC *bc, const CINTEnvVars *envs)
1696 {
1697 CINTg0_2e_2d(g, bc, envs);
1698 CINTg0_kj2d_4d(g, envs);
1699 }
CINTg0_2e_ik2d4d(double * g,struct _BC * bc,const CINTEnvVars * envs)1700 void CINTg0_2e_ik2d4d(double *g, struct _BC *bc, const CINTEnvVars *envs)
1701 {
1702 CINTg0_2e_2d(g, bc, envs);
1703 CINTg0_ik2d_4d(g, envs);
1704 }
CINTg0_2e_il2d4d(double * g,struct _BC * bc,const CINTEnvVars * envs)1705 void CINTg0_2e_il2d4d(double *g, struct _BC *bc, const CINTEnvVars *envs)
1706 {
1707 CINTg0_2e_2d(g, bc, envs);
1708 CINTg0_il2d_4d(g, envs);
1709 }
1710
1711 /*
1712 * g[i,k,l,j] = < ik | lj > = ( i j | k l )
1713 */
CINTg0_2e(double * g,const CINTEnvVars * envs)1714 FINT CINTg0_2e(double *g, const CINTEnvVars *envs)
1715 {
1716 FINT irys;
1717 FINT nroots = envs->nrys_roots;
1718 double aij = envs->ai[0] + envs->aj[0];
1719 double akl = envs->ak[0] + envs->al[0];
1720 double a0, a1, fac1, x;
1721 double u[MXRYSROOTS];
1722 double *w = g + envs->g_size * 2; // ~ gz
1723 double rijrkl[3];
1724 rijrkl[0] = envs->rij[0] - envs->rkl[0];
1725 rijrkl[1] = envs->rij[1] - envs->rkl[1];
1726 rijrkl[2] = envs->rij[2] - envs->rkl[2];
1727
1728 a1 = aij * akl;
1729 a0 = a1 / (aij + akl);
1730
1731 #ifdef WITH_RANGE_COULOMB
1732 const double omega = envs->env[PTR_RANGE_OMEGA];
1733 double theta = 0;
1734 if (omega != 0) {
1735 theta = omega * omega / (omega * omega + a0);
1736 if (omega > 0) { // long-range part of range-separated Coulomb
1737 a0 *= theta;
1738 }
1739 }
1740 #endif
1741
1742 x = a0 *(rijrkl[0] * rijrkl[0]
1743 + rijrkl[1] * rijrkl[1]
1744 + rijrkl[2] * rijrkl[2]);
1745
1746 #ifdef WITH_RANGE_COULOMB
1747 if (omega < 0) { // short-range part of range-separated Coulomb
1748 // FIXME:
1749 // very small erfc() leads to ~0 weights. They can cause
1750 // numerical issue in sr_rys_roots Use this cutoff as a
1751 // temporary solution to avoid the numerical issue
1752 double temp_cutoff = MIN(envs->expcutoff, EXPCUTOFF_SR - nroots);
1753 if (theta * x > temp_cutoff) {
1754 return 0;
1755 }
1756 CINTsr_rys_roots(nroots, x, sqrt(theta), u, w);
1757 } else {
1758 CINTrys_roots(nroots, x, u, w);
1759 if (omega > 0) {
1760 /* u[:] = tau^2 / (1 - tau^2)
1761 * omega^2u^2 = a0 * tau^2 / (theta^-1 - tau^2)
1762 * transform u[:] to theta^-1 tau^2 / (theta^-1 - tau^2)
1763 * so the rest code can be reused.
1764 */
1765 for (irys = 0; irys < nroots; irys++) {
1766 u[irys] /= u[irys] + 1 - u[irys] * theta;
1767 }
1768 }
1769 }
1770 #else
1771 CINTrys_roots(nroots, x, u, w);
1772 #endif
1773 fac1 = sqrt(a0 / (a1 * a1 * a1)) * envs->fac[0];
1774 if (envs->g_size == 1) {
1775 g[0] = 1;
1776 g[1] = 1;
1777 g[2] *= fac1;
1778 return 1;
1779 }
1780
1781 double u2, tmp1, tmp2, tmp3, tmp4, tmp5;
1782 double rijrx[3];
1783 double rklrx[3];
1784 rijrx[0] = envs->rij[0] - envs->rx_in_rijrx[0];
1785 rijrx[1] = envs->rij[1] - envs->rx_in_rijrx[1];
1786 rijrx[2] = envs->rij[2] - envs->rx_in_rijrx[2];
1787 rklrx[0] = envs->rkl[0] - envs->rx_in_rklrx[0];
1788 rklrx[1] = envs->rkl[1] - envs->rx_in_rklrx[1];
1789 rklrx[2] = envs->rkl[2] - envs->rx_in_rklrx[2];
1790 struct _BC bc;
1791 double *c00 = bc.c00;
1792 double *c0p = bc.c0p;
1793 double *b00 = bc.b00;
1794 double *b10 = bc.b10;
1795 double *b01 = bc.b01;
1796
1797 for (irys = 0; irys < nroots; irys++, c00+=3, c0p+=3) {
1798 /*
1799 *u(irys) = t2/(1-t2)
1800 *t2 = u(irys)/(1+u(irys))
1801 *u2 = aij*akl/(aij+akl)*t2/(1-t2)
1802 */
1803 u2 = a0 * u[irys];
1804 tmp4 = .5 / (u2 * (aij + akl) + a1);
1805 tmp5 = u2 * tmp4;
1806 tmp1 = 2. * tmp5;
1807 tmp2 = tmp1 * akl;
1808 tmp3 = tmp1 * aij;
1809 b00[irys] = tmp5;
1810 b10[irys] = tmp5 + tmp4 * akl;
1811 b01[irys] = tmp5 + tmp4 * aij;
1812 c00[0] = rijrx[0] - tmp2 * rijrkl[0];
1813 c00[1] = rijrx[1] - tmp2 * rijrkl[1];
1814 c00[2] = rijrx[2] - tmp2 * rijrkl[2];
1815 c0p[0] = rklrx[0] + tmp3 * rijrkl[0];
1816 c0p[1] = rklrx[1] + tmp3 * rijrkl[1];
1817 c0p[2] = rklrx[2] + tmp3 * rijrkl[2];
1818 w[irys] *= fac1;
1819 }
1820
1821 (*envs->f_g0_2d4d)(g, &bc, envs);
1822
1823 return 1;
1824 }
1825
1826 /*
1827 * ( \nabla i j | kl )
1828 */
CINTnabla1i_2e(double * f,const double * g,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)1829 void CINTnabla1i_2e(double *f, const double *g,
1830 const FINT li, const FINT lj, const FINT lk, const FINT ll,
1831 const CINTEnvVars *envs)
1832 {
1833 FINT i, j, k, l, n, ptr;
1834 const FINT di = envs->g_stride_i;
1835 const FINT dk = envs->g_stride_k;
1836 const FINT dl = envs->g_stride_l;
1837 const FINT dj = envs->g_stride_j;
1838 const FINT nroots = envs->nrys_roots;
1839 const double ai2 = -2 * envs->ai[0];
1840 DEF_GXYZ(const double, g, gx, gy, gz);
1841 DEF_GXYZ(double, f, fx, fy, fz);
1842
1843 const double *p1x = gx - di;
1844 const double *p1y = gy - di;
1845 const double *p1z = gz - di;
1846 const double *p2x = gx + di;
1847 const double *p2y = gy + di;
1848 const double *p2z = gz + di;
1849 for (j = 0; j <= lj; j++)
1850 for (l = 0; l <= ll; l++)
1851 for (k = 0; k <= lk; k++) {
1852 ptr = dj * j + dl * l + dk * k;
1853 //f(...,0,...) = -2*ai*g(...,1,...)
1854 for (n = ptr; n < ptr+nroots; n++) {
1855 fx[n] = ai2 * p2x[n];
1856 fy[n] = ai2 * p2y[n];
1857 fz[n] = ai2 * p2z[n];
1858 }
1859 ptr += di;
1860 //f(...,i,...) = i*g(...,i-1,...)-2*ai*g(...,i+1,...)
1861 for (i = 1; i <= li; i++) {
1862 for (n = ptr; n < ptr+nroots; n++) {
1863 fx[n] = i*p1x[n] + ai2*p2x[n];
1864 fy[n] = i*p1y[n] + ai2*p2y[n];
1865 fz[n] = i*p1z[n] + ai2*p2z[n];
1866 }
1867 ptr += di;
1868 }
1869 }
1870 }
1871
1872
1873 /*
1874 * ( i \nabla j | kl )
1875 */
CINTnabla1j_2e(double * f,const double * g,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)1876 void CINTnabla1j_2e(double *f, const double *g,
1877 const FINT li, const FINT lj, const FINT lk, const FINT ll,
1878 const CINTEnvVars *envs)
1879 {
1880 FINT i, j, k, l, n, ptr;
1881 const FINT di = envs->g_stride_i;
1882 const FINT dk = envs->g_stride_k;
1883 const FINT dl = envs->g_stride_l;
1884 const FINT dj = envs->g_stride_j;
1885 const FINT nroots = envs->nrys_roots;
1886 const double aj2 = -2 * envs->aj[0];
1887 DEF_GXYZ(const double, g, gx, gy, gz);
1888 DEF_GXYZ(double, f, fx, fy, fz);
1889
1890 const double *p1x = gx - dj;
1891 const double *p1y = gy - dj;
1892 const double *p1z = gz - dj;
1893 const double *p2x = gx + dj;
1894 const double *p2y = gy + dj;
1895 const double *p2z = gz + dj;
1896 //f(...,0,...) = -2*aj*g(...,1,...)
1897 for (l = 0; l <= ll; l++) {
1898 for (k = 0; k <= lk; k++) {
1899 ptr = dl * l + dk * k;
1900 for (i = 0; i <= li; i++) {
1901 for (n = ptr; n < ptr+nroots; n++) {
1902 fx[n] = aj2 * p2x[n];
1903 fy[n] = aj2 * p2y[n];
1904 fz[n] = aj2 * p2z[n];
1905 }
1906 ptr += di;
1907 }
1908 } }
1909 //f(...,j,...) = j*g(...,j-1,...)-2*aj*g(...,j+1,...)
1910 for (j = 1; j <= lj; j++) {
1911 for (l = 0; l <= ll; l++) {
1912 for (k = 0; k <= lk; k++) {
1913 ptr = dj * j + dl * l + dk * k;
1914 for (i = 0; i <= li; i++) {
1915 for (n = ptr; n < ptr+nroots; n++) {
1916 fx[n] = j*p1x[n] + aj2*p2x[n];
1917 fy[n] = j*p1y[n] + aj2*p2y[n];
1918 fz[n] = j*p1z[n] + aj2*p2z[n];
1919 }
1920 ptr += di;
1921 }
1922 } }
1923 }
1924 }
1925
1926
1927 /*
1928 * ( ij | \nabla k l )
1929 */
CINTnabla1k_2e(double * f,const double * g,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)1930 void CINTnabla1k_2e(double *f, const double *g,
1931 const FINT li, const FINT lj, const FINT lk, const FINT ll,
1932 const CINTEnvVars *envs)
1933 {
1934 FINT i, j, k, l, n, ptr;
1935 const FINT di = envs->g_stride_i;
1936 const FINT dk = envs->g_stride_k;
1937 const FINT dl = envs->g_stride_l;
1938 const FINT dj = envs->g_stride_j;
1939 const FINT nroots = envs->nrys_roots;
1940 const double ak2 = -2 * envs->ak[0];
1941 DEF_GXYZ(const double, g, gx, gy, gz);
1942 DEF_GXYZ(double, f, fx, fy, fz);
1943
1944 const double *p1x = gx - dk;
1945 const double *p1y = gy - dk;
1946 const double *p1z = gz - dk;
1947 const double *p2x = gx + dk;
1948 const double *p2y = gy + dk;
1949 const double *p2z = gz + dk;
1950 for (j = 0; j <= lj; j++)
1951 for (l = 0; l <= ll; l++) {
1952 ptr = dj * j + dl * l;
1953 //f(...,0,...) = -2*ak*g(...,1,...)
1954 for (i = 0; i <= li; i++) {
1955 for (n = ptr; n < ptr+nroots; n++) {
1956 fx[n] = ak2 * p2x[n];
1957 fy[n] = ak2 * p2y[n];
1958 fz[n] = ak2 * p2z[n];
1959 }
1960 ptr += di;
1961 }
1962 //f(...,k,...) = k*g(...,k-1,...)-2*ak*g(...,k+1,...)
1963 for (k = 1; k <= lk; k++) {
1964 ptr = dj * j + dl * l + dk * k;
1965 for (i = 0; i <= li; i++) {
1966 for (n = ptr; n < ptr+nroots; n++) {
1967 fx[n] = k*p1x[n] + ak2*p2x[n];
1968 fy[n] = k*p1y[n] + ak2*p2y[n];
1969 fz[n] = k*p1z[n] + ak2*p2z[n];
1970 }
1971 ptr += di;
1972 }
1973 }
1974 }
1975 }
1976
1977
1978 /*
1979 * ( ij | k \nabla l )
1980 */
CINTnabla1l_2e(double * f,const double * g,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)1981 void CINTnabla1l_2e(double *f, const double *g,
1982 const FINT li, const FINT lj, const FINT lk, const FINT ll,
1983 const CINTEnvVars *envs)
1984 {
1985 FINT i, j, k, l, n, ptr;
1986 const FINT di = envs->g_stride_i;
1987 const FINT dk = envs->g_stride_k;
1988 const FINT dl = envs->g_stride_l;
1989 const FINT dj = envs->g_stride_j;
1990 const FINT nroots = envs->nrys_roots;
1991 const double al2 = -2 * envs->al[0];
1992 DEF_GXYZ(const double, g, gx, gy, gz);
1993 DEF_GXYZ(double, f, fx, fy, fz);
1994
1995 const double *p1x = gx - dl;
1996 const double *p1y = gy - dl;
1997 const double *p1z = gz - dl;
1998 const double *p2x = gx + dl;
1999 const double *p2y = gy + dl;
2000 const double *p2z = gz + dl;
2001 for (j = 0; j <= lj; j++) {
2002 //f(...,0,...) = -2*al*g(...,1,...)
2003 for (k = 0; k <= lk; k++) {
2004 ptr = dj * j + dk * k;
2005 for (i = 0; i <= li; i++) {
2006 for (n = ptr; n < ptr+nroots; n++) {
2007 fx[n] = al2 * p2x[n];
2008 fy[n] = al2 * p2y[n];
2009 fz[n] = al2 * p2z[n];
2010 }
2011 ptr += di;
2012 }
2013 }
2014 //f(...,l,...) = l*g(...,l-1,...)-2*al*g(...,l+1,...)
2015 for (l = 1; l <= ll; l++) {
2016 for (k = 0; k <= lk; k++) {
2017 ptr = dj * j + dl * l + dk * k;
2018 for (i = 0; i <= li; i++, ptr += di) {
2019 for (n = ptr; n < ptr+nroots; n++) {
2020 fx[n] = l*p1x[n] + al2*p2x[n];
2021 fy[n] = l*p1y[n] + al2*p2y[n];
2022 fz[n] = l*p1z[n] + al2*p2z[n];
2023 } }
2024 }
2025 }
2026 }
2027 }
2028
2029 /*
2030 * ( x^1 i j | kl )
2031 * ri is the shift from the center R_O to the center of |i>
2032 * r - R_O = (r-R_i) + ri, ri = R_i - R_O
2033 */
CINTx1i_2e(double * f,const double * g,const double * ri,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)2034 void CINTx1i_2e(double *f, const double *g, const double *ri,
2035 const FINT li, const FINT lj, const FINT lk, const FINT ll,
2036 const CINTEnvVars *envs)
2037 {
2038 FINT i, j, k, l, n, ptr;
2039 const FINT di = envs->g_stride_i;
2040 const FINT dk = envs->g_stride_k;
2041 const FINT dl = envs->g_stride_l;
2042 const FINT dj = envs->g_stride_j;
2043 const FINT nroots = envs->nrys_roots;
2044 DEF_GXYZ(const double, g, gx, gy, gz);
2045 DEF_GXYZ(double, f, fx, fy, fz);
2046
2047 const double *p1x = gx + di;
2048 const double *p1y = gy + di;
2049 const double *p1z = gz + di;
2050 for (j = 0; j <= lj; j++)
2051 for (l = 0; l <= ll; l++) {
2052 for (k = 0; k <= lk; k++) {
2053 //f(...,0:li,...) = g(...,1:li+1,...) + ri(1)*g(...,0:li,...)
2054 ptr = dj * j + dl * l + dk * k;
2055 for (i = 0; i <= li; i++) {
2056 for (n = ptr; n < ptr+nroots; n++) {
2057 fx[n] = p1x[n] + ri[0] * gx[n];
2058 fy[n] = p1y[n] + ri[1] * gy[n];
2059 fz[n] = p1z[n] + ri[2] * gz[n];
2060 }
2061 ptr += di;
2062 }
2063 } }
2064 }
2065
2066
2067 /*
2068 * ( i x^1 j | kl )
2069 */
CINTx1j_2e(double * f,const double * g,const double * rj,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)2070 void CINTx1j_2e(double *f, const double *g, const double *rj,
2071 const FINT li, const FINT lj, const FINT lk, const FINT ll,
2072 const CINTEnvVars *envs)
2073 {
2074 FINT i, j, k, l, n, ptr;
2075 const FINT di = envs->g_stride_i;
2076 const FINT dk = envs->g_stride_k;
2077 const FINT dl = envs->g_stride_l;
2078 const FINT dj = envs->g_stride_j;
2079 const FINT nroots = envs->nrys_roots;
2080 DEF_GXYZ(const double, g, gx, gy, gz);
2081 DEF_GXYZ(double, f, fx, fy, fz);
2082
2083 const double *p1x = gx + dj;
2084 const double *p1y = gy + dj;
2085 const double *p1z = gz + dj;
2086 for (j = 0; j <= lj; j++)
2087 for (l = 0; l <= ll; l++) {
2088 for (k = 0; k <= lk; k++) {
2089 // f(...,0:lj,...) = g(...,1:lj+1,...) + rj(1)*g(...,0:lj,...)
2090 ptr = dj * j + dl * l + dk * k;
2091 for (i = 0; i <= li; i++) {
2092 for (n = ptr; n < ptr+nroots; n++) {
2093 fx[n] = p1x[n] + rj[0] * gx[n];
2094 fy[n] = p1y[n] + rj[1] * gy[n];
2095 fz[n] = p1z[n] + rj[2] * gz[n];
2096 }
2097 ptr += di;
2098 }
2099 } }
2100 }
2101
2102
2103 /*
2104 * ( ij | x^1 k l )
2105 */
CINTx1k_2e(double * f,const double * g,const double * rk,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)2106 void CINTx1k_2e(double *f, const double *g, const double *rk,
2107 const FINT li, const FINT lj, const FINT lk, const FINT ll,
2108 const CINTEnvVars *envs)
2109 {
2110 FINT i, j, k, l, n, ptr;
2111 const FINT di = envs->g_stride_i;
2112 const FINT dk = envs->g_stride_k;
2113 const FINT dl = envs->g_stride_l;
2114 const FINT dj = envs->g_stride_j;
2115 const FINT nroots = envs->nrys_roots;
2116 DEF_GXYZ(const double, g, gx, gy, gz);
2117 DEF_GXYZ(double, f, fx, fy, fz);
2118
2119 const double *p1x = gx + dk;
2120 const double *p1y = gy + dk;
2121 const double *p1z = gz + dk;
2122 for (j = 0; j <= lj; j++)
2123 for (l = 0; l <= ll; l++) {
2124 for (k = 0; k <= lk; k++) {
2125 // f(...,0:lk,...) = g(...,1:lk+1,...) + rk(1)*g(...,0:lk,...)
2126 ptr = dj * j + dl * l + dk * k;
2127 for (i = 0; i <= li; i++) {
2128 for (n = ptr; n < ptr+nroots; n++) {
2129 fx[n] = p1x[n] + rk[0] * gx[n];
2130 fy[n] = p1y[n] + rk[1] * gy[n];
2131 fz[n] = p1z[n] + rk[2] * gz[n];
2132 }
2133 ptr += di;
2134 }
2135 } }
2136 }
2137
2138
2139 /*
2140 * ( i j | x^1 kl )
2141 */
CINTx1l_2e(double * f,const double * g,const double * rl,const FINT li,const FINT lj,const FINT lk,const FINT ll,const CINTEnvVars * envs)2142 void CINTx1l_2e(double *f, const double *g, const double *rl,
2143 const FINT li, const FINT lj, const FINT lk, const FINT ll,
2144 const CINTEnvVars *envs)
2145 {
2146 FINT i, j, k, l, n, ptr;
2147 const FINT di = envs->g_stride_i;
2148 const FINT dk = envs->g_stride_k;
2149 const FINT dl = envs->g_stride_l;
2150 const FINT dj = envs->g_stride_j;
2151 const FINT nroots = envs->nrys_roots;
2152 DEF_GXYZ(const double, g, gx, gy, gz);
2153 DEF_GXYZ(double, f, fx, fy, fz);
2154
2155 const double *p1x = gx + dl;
2156 const double *p1y = gy + dl;
2157 const double *p1z = gz + dl;
2158 for (j = 0; j <= lj; j++)
2159 for (l = 0; l <= ll; l++) {
2160 for (k = 0; k <= lk; k++) {
2161 // f(...,0:ll,...) = g(...,1:ll+1,...) + rl(1)*g(...,0:ll,...)
2162 ptr = dj * j + dl * l + dk * k;
2163 for (i = 0; i <= li; i++) {
2164 for (n = ptr; n < ptr+nroots; n++) {
2165 fx[n] = p1x[n] + rl[0] * gx[n];
2166 fy[n] = p1y[n] + rl[1] * gy[n];
2167 fz[n] = p1z[n] + rl[2] * gz[n];
2168 }
2169 ptr += di;
2170 }
2171 } }
2172 }
2173
2174