1 /*
2 * Copyright (C) 2013- Qiming Sun <osirpt.sun@gmail.com>
3 *
4 * 3-center 2-electron integrals
5 */
6
7 #include <stdlib.h>
8 #include <string.h>
9 #include <math.h>
10 #include "cint_bas.h"
11 #include "optimizer.h"
12 #include "g2e.h"
13 #include "cint2e.h"
14 #include "misc.h"
15 #include "cart2sph.h"
16 #include "c2f.h"
17
18 #define gctrg gout
19 #define gctrm gctr
20 #define mempty empty
21 #define m_ctr n_comp
22 #define ALIAS_ADDR_IF_EQUAL(x, y) \
23 if (y##_ctr == 1) { \
24 gctr##x = gctr##y; \
25 x##empty = y##empty; \
26 } else { \
27 gctr##x = g1; \
28 g1 += len##x; \
29 }
30
31 #define PRIM2CTR0(ctrsymb, gp, ngp) \
32 if (ctrsymb##_ctr > 1) {\
33 if (*ctrsymb##empty) { \
34 CINTprim_to_ctr_0(gctr##ctrsymb, gp, c##ctrsymb+ctrsymb##p, \
35 ngp, ctrsymb##_prim, ctrsymb##_ctr, \
36 non0ctr##ctrsymb[ctrsymb##p], \
37 non0idx##ctrsymb+ctrsymb##p*ctrsymb##_ctr); \
38 } else { \
39 CINTprim_to_ctr_1(gctr##ctrsymb, gp, c##ctrsymb+ctrsymb##p, \
40 ngp, ctrsymb##_prim, ctrsymb##_ctr, \
41 non0ctr##ctrsymb[ctrsymb##p], \
42 non0idx##ctrsymb+ctrsymb##p*ctrsymb##_ctr); \
43 } \
44 } \
45 *ctrsymb##empty = 0
46
47 #define TRANSPOSE(a) \
48 if (*empty) { \
49 CINTdmat_transpose(gctr, a, nf*nc, n_comp); \
50 } else { \
51 CINTdplus_transpose(gctr, a, nf*nc, n_comp); \
52 } \
53 *empty = 0;
54
55
CINT3c2e_loop_nopt(double * gctr,CINTEnvVars * envs,double * cache,FINT * empty)56 FINT CINT3c2e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache, FINT *empty)
57 {
58 FINT *shls = envs->shls;
59 FINT *bas = envs->bas;
60 double *env = envs->env;
61 FINT i_sh = shls[0];
62 FINT j_sh = shls[1];
63 FINT k_sh = shls[2];
64 FINT i_ctr = envs->x_ctr[0];
65 FINT j_ctr = envs->x_ctr[1];
66 FINT k_ctr = envs->x_ctr[2];
67 FINT i_prim = bas(NPRIM_OF, i_sh);
68 FINT j_prim = bas(NPRIM_OF, j_sh);
69 FINT k_prim = bas(NPRIM_OF, k_sh);
70 //double *ri = envs->ri;
71 //double *rj = envs->rj;
72 double *ai = env + bas(PTR_EXP, i_sh);
73 double *aj = env + bas(PTR_EXP, j_sh);
74 double *ak = env + bas(PTR_EXP, k_sh);
75 double *ci = env + bas(PTR_COEFF, i_sh);
76 double *cj = env + bas(PTR_COEFF, j_sh);
77 double *ck = env + bas(PTR_COEFF, k_sh);
78
79 double expcutoff = envs->expcutoff;
80 double *log_maxci, *log_maxcj;
81 PairData *pdata_base, *pdata_ij;
82 MALLOC_INSTACK(log_maxci, i_prim+j_prim);
83 MALLOC_INSTACK(pdata_base, i_prim*j_prim);
84 log_maxcj = log_maxci + i_prim;
85 CINTOpt_log_max_pgto_coeff(log_maxci, ci, i_prim, i_ctr);
86 CINTOpt_log_max_pgto_coeff(log_maxcj, cj, j_prim, j_ctr);
87 if (CINTset_pairdata(pdata_base, ai, aj, envs->ri, envs->rj,
88 log_maxci, log_maxcj, envs->li_ceil, envs->lj_ceil,
89 i_prim, j_prim, SQUARE(envs->rirj), expcutoff)) {
90 return 0;
91 }
92
93 FINT n_comp = envs->ncomp_e1 * envs->ncomp_tensor;
94 size_t nf = envs->nf;
95 double fac1i, fac1j, fac1k;
96 FINT ip, jp, kp;
97 FINT _empty[4] = {1, 1, 1, 1};
98 FINT *iempty = _empty + 0;
99 FINT *jempty = _empty + 1;
100 FINT *kempty = _empty + 2;
101 FINT *gempty = _empty + 3;
102 /* COMMON_ENVS_AND_DECLARE end */
103
104 double expij;
105 double *rij;
106 //double dist_ij = SQUARE(envs->rirj);
107
108 FINT *idx;
109 MALLOC_INSTACK(idx, nf * 3);
110 CINTg2e_index_xyz(idx, envs);
111
112 FINT *non0ctri, *non0ctrj, *non0ctrk;
113 FINT *non0idxi, *non0idxj, *non0idxk;
114 MALLOC_INSTACK(non0ctri, i_prim+j_prim+k_prim+i_prim*i_ctr+j_prim*j_ctr+k_prim*k_ctr);
115 non0ctrj = non0ctri + i_prim;
116 non0ctrk = non0ctrj + j_prim;
117 non0idxi = non0ctrk + k_prim;
118 non0idxj = non0idxi + i_prim*i_ctr;
119 non0idxk = non0idxj + j_prim*j_ctr;
120 CINTOpt_non0coeff_byshell(non0idxi, non0ctri, ci, i_prim, i_ctr);
121 CINTOpt_non0coeff_byshell(non0idxj, non0ctrj, cj, j_prim, j_ctr);
122 CINTOpt_non0coeff_byshell(non0idxk, non0ctrk, ck, k_prim, k_ctr);
123
124 FINT nc = i_ctr * j_ctr * k_ctr;
125 // (irys,i,j,k,l,coord,0:1); +1 for nabla-r12
126 size_t leng = envs->g_size * 3 * ((1<<envs->gbits)+1);
127 size_t lenk = nf * nc * n_comp; // gctrk
128 size_t lenj = nf * i_ctr * j_ctr * n_comp; // gctrj
129 size_t leni = nf * i_ctr * n_comp; // gctri
130 size_t len0 = nf * n_comp; // gout
131 size_t len = leng + lenk + lenj + leni + len0;
132 double *g;
133 MALLOC_INSTACK(g, len); // must be allocated last in this function
134 double *g1 = g + leng;
135 double *gout, *gctri, *gctrj, *gctrk;
136
137 ALIAS_ADDR_IF_EQUAL(k, m);
138 ALIAS_ADDR_IF_EQUAL(j, k);
139 ALIAS_ADDR_IF_EQUAL(i, j);
140 ALIAS_ADDR_IF_EQUAL(g, i);
141
142 for (kp = 0; kp < k_prim; kp++) {
143 envs->ak[0] = ak[kp];
144 if (k_ctr == 1) {
145 fac1k = envs->common_factor * ck[kp];
146 } else {
147 fac1k = envs->common_factor;
148 *jempty = 1;
149 }
150
151 pdata_ij = pdata_base;
152 for (jp = 0; jp < j_prim; jp++) {
153 envs->aj[0] = aj[jp];
154 if (j_ctr == 1) {
155 fac1j = fac1k * cj[jp];
156 } else {
157 fac1j = fac1k;
158 *iempty = 1;
159 }
160 for (ip = 0; ip < i_prim; ip++, pdata_ij++) {
161 if (pdata_ij->cceij > expcutoff) {
162 goto i_contracted;
163 }
164 envs->ai[0] = ai[ip];
165 expij = pdata_ij->eij;
166 rij = pdata_ij->rij;
167 envs->rij[0] = rij[0];
168 envs->rij[1] = rij[1];
169 envs->rij[2] = rij[2];
170 if (i_ctr == 1) {
171 fac1i = fac1j*ci[ip]*expij;
172 } else {
173 fac1i = fac1j*expij;
174 }
175 envs->fac[0] = fac1i;
176 if ((*envs->f_g0_2e)(g, envs)) {
177 (*envs->f_gout)(gout, g, idx, envs, *gempty);
178 PRIM2CTR0(i, gout, len0);
179 }
180 i_contracted: ;
181 } // end loop i_prim
182 if (!*iempty) {
183 PRIM2CTR0(j, gctri, leni);
184 }
185 } // end loop j_prim
186 if (!*jempty) {
187 PRIM2CTR0(k, gctrj, lenj);
188 }
189 } // end loop k_prim
190
191 if (n_comp > 1 && !*kempty) {
192 TRANSPOSE(gctrk);
193 }
194 return !*empty;
195 }
196
197
198 #define COMMON_ENVS_AND_DECLARE \
199 FINT *shls = envs->shls; \
200 FINT *bas = envs->bas; \
201 double *env = envs->env; \
202 FINT i_sh = shls[0]; \
203 FINT j_sh = shls[1]; \
204 CINTOpt *opt = envs->opt; \
205 if (opt->pairdata != NULL && \
206 opt->pairdata[i_sh*opt->nbas+j_sh] == NOVALUE) { \
207 return 0; \
208 } \
209 FINT k_sh = shls[2]; \
210 FINT i_ctr = envs->x_ctr[0]; \
211 FINT j_ctr = envs->x_ctr[1]; \
212 FINT k_ctr = envs->x_ctr[2]; \
213 FINT i_prim = bas(NPRIM_OF, i_sh); \
214 FINT j_prim = bas(NPRIM_OF, j_sh); \
215 FINT k_prim = bas(NPRIM_OF, k_sh); \
216 double *ai = env + bas(PTR_EXP, i_sh); \
217 double *aj = env + bas(PTR_EXP, j_sh); \
218 double *ak = env + bas(PTR_EXP, k_sh); \
219 double *ci = env + bas(PTR_COEFF, i_sh); \
220 double *cj = env + bas(PTR_COEFF, j_sh); \
221 double *ck = env + bas(PTR_COEFF, k_sh); \
222 double expcutoff = envs->expcutoff; \
223 PairData *pdata_base, *pdata_ij; \
224 if (opt->pairdata != NULL) { \
225 pdata_base = opt->pairdata[i_sh*opt->nbas+j_sh]; \
226 } else { \
227 double *log_maxci = opt->log_max_coeff[i_sh]; \
228 double *log_maxcj = opt->log_max_coeff[j_sh]; \
229 MALLOC_INSTACK(pdata_base, i_prim*j_prim); \
230 if (CINTset_pairdata(pdata_base, ai, aj, envs->ri, envs->rj, \
231 log_maxci, log_maxcj, envs->li_ceil, envs->lj_ceil, \
232 i_prim, j_prim, SQUARE(envs->rirj), expcutoff)) { \
233 return 0; \
234 } \
235 } \
236 FINT n_comp = envs->ncomp_e1 * envs->ncomp_tensor; \
237 size_t nf = envs->nf; \
238 double fac1i, fac1j, fac1k; \
239 FINT ip, jp, kp; \
240 FINT _empty[4] = {1, 1, 1, 1}; \
241 FINT *iempty = _empty + 0; \
242 FINT *jempty = _empty + 1; \
243 FINT *kempty = _empty + 2; \
244 FINT *gempty = _empty + 3; \
245 FINT *non0ctri = opt->non0ctr[i_sh]; \
246 FINT *non0ctrj = opt->non0ctr[j_sh]; \
247 FINT *non0idxi = opt->sortedidx[i_sh]; \
248 FINT *non0idxj = opt->sortedidx[j_sh]; \
249 FINT *non0ctrk, *non0idxk; \
250 MALLOC_INSTACK(non0ctrk, k_prim+k_prim*k_ctr); \
251 non0idxk = non0ctrk + k_prim; \
252 CINTOpt_non0coeff_byshell(non0idxk, non0ctrk, ck, k_prim, k_ctr); \
253 double expij; \
254 double *rij; \
255 FINT *idx = opt->index_xyz_array[envs->i_l*LMAX1*LMAX1 \
256 +envs->j_l*LMAX1+envs->k_l]; \
257 if (idx == NULL) { \
258 MALLOC_INSTACK(idx, nf * 3); \
259 CINTg2e_index_xyz(idx, envs); \
260 }
261
262 #define SET_RIJ \
263 if (pdata_ij->cceij > expcutoff) { \
264 goto i_contracted; \
265 } \
266 envs->ai[0] = ai[ip]; \
267 expij = pdata_ij->eij; \
268 rij = pdata_ij->rij; \
269 envs->rij[0] = rij[0]; \
270 envs->rij[1] = rij[1]; \
271 envs->rij[2] = rij[2];
272
273 #define PRIM2CTR(ctrsymb, gp, ngp) \
274 if (ctrsymb##_ctr > 1) {\
275 if (*ctrsymb##empty) { \
276 CINTprim_to_ctr_0(gctr##ctrsymb, gp, c##ctrsymb+ctrsymb##p, \
277 ngp, ctrsymb##_prim, ctrsymb##_ctr, \
278 non0ctr##ctrsymb[ctrsymb##p], \
279 non0idx##ctrsymb+ctrsymb##p*ctrsymb##_ctr); \
280 } else { \
281 CINTprim_to_ctr_1(gctr##ctrsymb, gp, c##ctrsymb+ctrsymb##p, \
282 ngp, ctrsymb##_prim, ctrsymb##_ctr, \
283 non0ctr##ctrsymb[ctrsymb##p], \
284 non0idx##ctrsymb+ctrsymb##p*ctrsymb##_ctr); \
285 } \
286 } \
287 *ctrsymb##empty = 0
288
289
290 // i_ctr = j_ctr = k_ctr = 1;
CINT3c2e_111_loop(double * gctr,CINTEnvVars * envs,double * cache,FINT * empty)291 FINT CINT3c2e_111_loop(double *gctr, CINTEnvVars *envs, double *cache, FINT *empty)
292 {
293 COMMON_ENVS_AND_DECLARE;
294 FINT nc = 1;
295 size_t leng = envs->g_size * 3 * ((1<<envs->gbits)+1);
296 size_t len0 = envs->nf * n_comp;
297 size_t len = leng + len0;
298 double *g;
299 MALLOC_INSTACK(g, len);
300 double *gout;
301 if (n_comp == 1) {
302 gout = gctr;
303 gempty = empty;
304 } else {
305 gout = g + leng;
306 }
307
308 for (kp = 0; kp < k_prim; kp++) {
309 envs->ak[0] = ak[kp];
310 fac1k = envs->common_factor * ck[kp];
311
312 pdata_ij = pdata_base;
313 for (jp = 0; jp < j_prim; jp++) {
314 envs->aj[0] = aj[jp];
315 fac1j = fac1k * cj[jp];
316 for (ip = 0; ip < i_prim; ip++, pdata_ij++) {
317 SET_RIJ;
318 fac1i = fac1j*ci[ip]*expij;
319 envs->fac[0] = fac1i;
320 if ((*envs->f_g0_2e)(g, envs)) {
321 (*envs->f_gout)(gout, g, idx, envs, *gempty);
322 *gempty = 0;
323 }
324 i_contracted: ;
325 } // end loop i_prim
326 } // end loop j_prim
327 } // end loop k_prim
328
329 if (n_comp > 1 && !*gempty) {
330 TRANSPOSE(gout);
331 }
332 return !*empty;
333 }
334
335 // i_ctr = n; j_ctr = k_ctr = 1;
CINT3c2e_n11_loop(double * gctr,CINTEnvVars * envs,double * cache,FINT * empty)336 FINT CINT3c2e_n11_loop(double *gctr, CINTEnvVars *envs, double *cache, FINT *empty)
337 {
338 COMMON_ENVS_AND_DECLARE;
339
340 FINT nc = i_ctr;
341 size_t leng = envs->g_size * 3 * ((1<<envs->gbits)+1);
342 size_t leni = nf * i_ctr * n_comp; // gctri
343 size_t len0 = nf * n_comp; // gout
344 size_t len = leng + leni + len0;
345 double *g;
346 MALLOC_INSTACK(g, len);
347 double *g1 = g + leng;
348 double *gout, *gctri;
349 ALIAS_ADDR_IF_EQUAL(i, m);
350 gout = g1;
351
352 for (kp = 0; kp < k_prim; kp++) {
353 envs->ak[0] = ak[kp];
354 fac1k = envs->common_factor * ck[kp];
355 pdata_ij = pdata_base;
356 for (jp = 0; jp < j_prim; jp++) {
357 envs->aj[0] = aj[jp];
358 fac1j = fac1k * cj[jp];
359 for (ip = 0; ip < i_prim; ip++, pdata_ij++) {
360 SET_RIJ;
361 fac1i = fac1j*expij;
362 envs->fac[0] = fac1i;
363 if ((*envs->f_g0_2e)(g, envs)) {
364 (*envs->f_gout)(gout, g, idx, envs, 1);
365 PRIM2CTR(i, gout, len0);
366 }
367 i_contracted: ;
368 } // end loop i_prim
369 } // end loop j_prim
370 } // end loop k_prim
371
372 if (n_comp > 1 && !*iempty) {
373 TRANSPOSE(gctri);
374 }
375 return !*empty;
376 }
377
378 // j_ctr = n; i_ctr = k_ctr = 1;
CINT3c2e_1n1_loop(double * gctr,CINTEnvVars * envs,double * cache,FINT * empty)379 FINT CINT3c2e_1n1_loop(double *gctr, CINTEnvVars *envs, double *cache, FINT *empty)
380 {
381 COMMON_ENVS_AND_DECLARE;
382
383 FINT nc = j_ctr;
384 size_t leng = envs->g_size * 3 * ((1<<envs->gbits)+1);
385 size_t lenj = nf * j_ctr * n_comp; // gctrj
386 size_t len0 = nf * n_comp; // gout
387 size_t len = leng + lenj + len0;
388 double *g;
389 MALLOC_INSTACK(g, len);
390 double *g1 = g + leng;
391 double *gout, *gctrj;
392 ALIAS_ADDR_IF_EQUAL(j, m);
393 gout = g1;
394
395 for (kp = 0; kp < k_prim; kp++) {
396 envs->ak[0] = ak[kp];
397 fac1k = envs->common_factor * ck[kp];
398 pdata_ij = pdata_base;
399 for (jp = 0; jp < j_prim; jp++) {
400 envs->aj[0] = aj[jp];
401 fac1j = fac1k;
402 *iempty = 1;
403 for (ip = 0; ip < i_prim; ip++, pdata_ij++) {
404 SET_RIJ;
405 fac1i = fac1j*ci[ip]*expij;
406 envs->fac[0] = fac1i;
407 if ((*envs->f_g0_2e)(g, envs)) {
408 (*envs->f_gout)(gout, g, idx, envs, *iempty);
409 *iempty = 0;
410 }
411 i_contracted: ;
412 } // end loop i_prim
413 if (!*iempty) {
414 PRIM2CTR(j, gout, len0);
415 }
416 } // end loop j_prim
417 } // end loop k_prim
418
419 if (n_comp > 1 && !*jempty) {
420 TRANSPOSE(gctrj);
421 }
422 return !*empty;
423 }
424
425
CINT3c2e_loop(double * gctr,CINTEnvVars * envs,double * cache,FINT * empty)426 FINT CINT3c2e_loop(double *gctr, CINTEnvVars *envs, double *cache, FINT *empty)
427 {
428 COMMON_ENVS_AND_DECLARE;
429 FINT nc = i_ctr * j_ctr * k_ctr;
430 // (irys,i,j,k,coord,0:1); +1 for nabla-r12
431 size_t leng = envs->g_size * 3 * ((1<<envs->gbits)+1);
432 size_t lenk = nf * nc * n_comp; // gctrk
433 size_t lenj = nf * i_ctr * j_ctr * n_comp; // gctrj
434 size_t leni = nf * i_ctr * n_comp; // gctri
435 size_t len0 = nf * n_comp; // gout
436 size_t len = leng + lenk + lenj + leni + len0;
437 double *g;
438 MALLOC_INSTACK(g, len);
439 double *g1 = g + leng;
440 double *gout, *gctri, *gctrj, *gctrk;
441
442 ALIAS_ADDR_IF_EQUAL(k, m);
443 ALIAS_ADDR_IF_EQUAL(j, k);
444 ALIAS_ADDR_IF_EQUAL(i, j);
445 ALIAS_ADDR_IF_EQUAL(g, i);
446
447 for (kp = 0; kp < k_prim; kp++) {
448 envs->ak[0] = ak[kp];
449 if (k_ctr == 1) {
450 fac1k = envs->common_factor * ck[kp];
451 } else {
452 fac1k = envs->common_factor;
453 *jempty = 1;
454 }
455 pdata_ij = pdata_base;
456 for (jp = 0; jp < j_prim; jp++) {
457 envs->aj[0] = aj[jp];
458 if (j_ctr == 1) {
459 fac1j = fac1k * cj[jp];
460 } else {
461 fac1j = fac1k;
462 *iempty = 1;
463 }
464 for (ip = 0; ip < i_prim; ip++, pdata_ij++) {
465 SET_RIJ;
466 if (i_ctr == 1) {
467 fac1i = fac1j*ci[ip]*expij;
468 } else {
469 fac1i = fac1j*expij;
470 }
471 envs->fac[0] = fac1i;
472 if ((*envs->f_g0_2e)(g, envs)) {
473 (*envs->f_gout)(gout, g, idx, envs, *gempty);
474 PRIM2CTR(i, gout, len0);
475 }
476 i_contracted: ;
477 } // end loop i_prim
478 if (!*iempty) {
479 PRIM2CTR(j, gctri, leni);
480 }
481 } // end loop j_prim
482 if (!*jempty) {
483 PRIM2CTR0(k, gctrj, lenj);
484 }
485 } // end loop k_prim
486
487 if (n_comp > 1 && !*kempty) {
488 TRANSPOSE(gctrk);
489 }
490 return !*empty;
491 }
492
493 static FINT (*CINTf_3c2e_loop[8])(double *, CINTEnvVars *, double *, FINT *) = {
494 CINT3c2e_loop,
495 CINT3c2e_loop,
496 CINT3c2e_loop,
497 CINT3c2e_n11_loop,
498 CINT3c2e_loop,
499 CINT3c2e_1n1_loop,
500 CINT3c2e_loop,
501 CINT3c2e_111_loop,
502 };
503
504 #define PAIRDATA_NON0IDX_SIZE(ps) \
505 FINT *bas = envs->bas; \
506 FINT *shls = envs->shls; \
507 FINT i_prim = bas(NPRIM_OF, shls[0]); \
508 FINT j_prim = bas(NPRIM_OF, shls[1]); \
509 FINT k_prim = bas(NPRIM_OF, shls[2]); \
510 FINT ps = (i_prim*j_prim * 5 \
511 + i_prim * x_ctr[0] \
512 + j_prim * x_ctr[1] \
513 + k_prim * x_ctr[2] \
514 +(i_prim+j_prim)*2 + k_prim + envs->nf*3 + 16);
515
CINT3c2e_drv(double * out,FINT * dims,CINTEnvVars * envs,CINTOpt * opt,double * cache,void (* f_e1_c2s)(),FINT is_ssc)516 CACHE_SIZE_T CINT3c2e_drv(double *out, FINT *dims, CINTEnvVars *envs, CINTOpt *opt,
517 double *cache, void (*f_e1_c2s)(), FINT is_ssc)
518 {
519 FINT *x_ctr = envs->x_ctr;
520 size_t nc = envs->nf * x_ctr[0] * x_ctr[1] * x_ctr[2];
521 FINT n_comp = envs->ncomp_e1 * envs->ncomp_tensor;
522 if (out == NULL) {
523 PAIRDATA_NON0IDX_SIZE(pdata_size);
524 CACHE_SIZE_T leng = envs->g_size*3*((1<<envs->gbits)+1);
525 CACHE_SIZE_T len0 = envs->nf*n_comp;
526 CACHE_SIZE_T cache_size = MAX(leng+len0+nc*n_comp*3 + pdata_size,
527 nc*n_comp+envs->nf*3);
528 return cache_size;
529 }
530 double *stack = NULL;
531 if (cache == NULL) {
532 PAIRDATA_NON0IDX_SIZE(pdata_size);
533 size_t leng = envs->g_size*3*((1<<envs->gbits)+1);
534 size_t len0 = envs->nf*n_comp;
535 size_t cache_size = MAX(leng+len0+nc*n_comp*3 + pdata_size,
536 nc*n_comp+envs->nf*3);
537 stack = malloc(sizeof(double)*cache_size);
538 cache = stack;
539 }
540 double *gctr;
541 MALLOC_INSTACK(gctr, nc*n_comp);
542
543 FINT n;
544 FINT empty = 1;
545 if (opt != NULL) {
546 envs->opt = opt;
547 n = ((envs->x_ctr[0]==1) << 2) + ((envs->x_ctr[1]==1) << 1) + (envs->x_ctr[2]==1);
548 CINTf_3c2e_loop[n](gctr, envs, cache, &empty);
549 } else {
550 CINT3c2e_loop_nopt(gctr, envs, cache, &empty);
551 }
552
553 FINT counts[4];
554 if (f_e1_c2s == &c2s_sph_3c2e1) {
555 counts[0] = (envs->i_l*2+1) * x_ctr[0];
556 counts[1] = (envs->j_l*2+1) * x_ctr[1];
557 if (is_ssc) {
558 counts[2] = envs->nfk * x_ctr[2];
559 } else {
560 counts[2] = (envs->k_l*2+1) * x_ctr[2];
561 }
562 } else {
563 counts[0] = envs->nfi * x_ctr[0];
564 counts[1] = envs->nfj * x_ctr[1];
565 counts[2] = envs->nfk * x_ctr[2];
566 }
567 counts[3] = 1;
568 if (dims == NULL) {
569 dims = counts;
570 }
571 FINT nout = dims[0] * dims[1] * dims[2];
572 if (!empty) {
573 for (n = 0; n < n_comp; n++) {
574 (*f_e1_c2s)(out+nout*n, gctr+nc*n, dims, envs, cache);
575 }
576 } else {
577 for (n = 0; n < n_comp; n++) {
578 c2s_dset0(out+nout*n, dims, counts);
579 }
580 }
581 if (stack != NULL) {
582 free(stack);
583 }
584 return !empty;
585 }
CINT3c2e_spinor_drv(double complex * out,FINT * dims,CINTEnvVars * envs,CINTOpt * opt,double * cache,void (* f_e1_c2s)(),FINT is_ssc)586 CACHE_SIZE_T CINT3c2e_spinor_drv(double complex *out, FINT *dims, CINTEnvVars *envs, CINTOpt *opt,
587 double *cache, void (*f_e1_c2s)(), FINT is_ssc)
588 {
589 FINT *x_ctr = envs->x_ctr;
590 FINT counts[4];
591 counts[0] = CINTcgto_spinor(envs->shls[0], envs->bas);
592 counts[1] = CINTcgto_spinor(envs->shls[1], envs->bas);
593 if (is_ssc) {
594 counts[2] = envs->nfk * x_ctr[2];
595 } else {
596 counts[2] = (envs->k_l*2+1) * x_ctr[2];
597 }
598 counts[3] = 1;
599 FINT nc = envs->nf * x_ctr[0] * x_ctr[1] * x_ctr[2];
600 FINT n_comp = envs->ncomp_e1 * envs->ncomp_e2 * envs->ncomp_tensor;
601 if (out == NULL) {
602 PAIRDATA_NON0IDX_SIZE(pdata_size);
603 size_t leng = envs->g_size*3*((1<<envs->gbits)+1);
604 size_t len0 = envs->nf*n_comp;
605 FINT cache_size = MAX(leng+len0+nc*n_comp*3 + pdata_size,
606 nc*n_comp + envs->nf*14*OF_CMPLX);
607 return cache_size;
608 }
609 double *stack = NULL;
610 if (cache == NULL) {
611 PAIRDATA_NON0IDX_SIZE(pdata_size);
612 size_t leng = envs->g_size*3*((1<<envs->gbits)+1);
613 size_t len0 = envs->nf*n_comp;
614 FINT cache_size = MAX(leng+len0+nc*n_comp*3 + pdata_size,
615 nc*n_comp + envs->nf*14*OF_CMPLX);
616 stack = malloc(sizeof(double)*cache_size);
617 cache = stack;
618 }
619 double *gctr;
620 MALLOC_INSTACK(gctr, nc*n_comp);
621
622 FINT n;
623 FINT empty = 1;
624 if (opt != NULL) {
625 envs->opt = opt;
626 n = ((envs->x_ctr[0]==1) << 2) + ((envs->x_ctr[1]==1) << 1) + (envs->x_ctr[2]==1);
627 CINTf_3c2e_loop[n](gctr, envs, cache, &empty);
628 } else {
629 CINT3c2e_loop_nopt(gctr, envs, cache, &empty);
630 }
631
632 if (dims == NULL) {
633 dims = counts;
634 }
635 FINT nout = dims[0] * dims[1] * dims[2];
636 if (!empty) {
637 for (n = 0; n < envs->ncomp_e2 * envs->ncomp_tensor; n++) {
638 (*f_e1_c2s)(out+nout*n, gctr, dims, envs, cache);
639 gctr += nc * envs->ncomp_e1;
640 }
641 } else {
642 for (n = 0; n < envs->ncomp_e2 * envs->ncomp_tensor; n++) {
643 c2s_zset0(out+nout*n, dims, counts);
644 }
645 }
646 if (stack != NULL) {
647 free(stack);
648 }
649 return !empty;
650 }
651
652
int3c2e_sph(double * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)653 CACHE_SIZE_T int3c2e_sph(double *out, FINT *dims, FINT *shls, FINT *atm, FINT natm,
654 FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache)
655 {
656 FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
657 CINTEnvVars envs;
658 CINTinit_int3c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
659 envs.f_gout = &CINTgout2e;
660 return CINT3c2e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c2e1, 0);
661 }
int3c2e_optimizer(CINTOpt ** opt,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env)662 void int3c2e_optimizer(CINTOpt **opt, FINT *atm, FINT natm,
663 FINT *bas, FINT nbas, double *env)
664 {
665 FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
666 CINTall_3c2e_optimizer(opt, ng, atm, natm, bas, nbas, env);
667 }
668
int3c2e_cart(double * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)669 CACHE_SIZE_T int3c2e_cart(double *out, FINT *dims, FINT *shls, FINT *atm, FINT natm,
670 FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache)
671 {
672 FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
673 CINTEnvVars envs;
674 CINTinit_int3c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
675 envs.f_gout = &CINTgout2e;
676 return CINT3c2e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c2e1, 0);
677 }
678
int3c2e_spinor(double complex * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)679 CACHE_SIZE_T int3c2e_spinor(double complex *out, FINT *dims, FINT *shls, FINT *atm, FINT natm,
680 FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache)
681 {
682 FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
683 CINTEnvVars envs;
684 CINTinit_int3c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
685 envs.f_gout = &CINTgout2e;
686 return CINT3c2e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 0);
687 }
688
int3c2e_sph_ssc(double * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)689 CACHE_SIZE_T int3c2e_sph_ssc(double *out, FINT *dims, FINT *shls, FINT *atm, FINT natm,
690 FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache)
691 {
692 FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
693 CINTEnvVars envs;
694 CINTinit_int3c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
695 envs.f_gout = &CINTgout2e;
696 return CINT3c2e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c2e1_ssc, 1);
697 }
int3c2e_spinor_ssc(double complex * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)698 CACHE_SIZE_T int3c2e_spinor_ssc(double complex *out, FINT *dims, FINT *shls, FINT *atm, FINT natm,
699 FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache)
700 {
701 FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
702 CINTEnvVars envs;
703 CINTinit_int3c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
704 envs.f_gout = &CINTgout2e;
705 return CINT3c2e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1_ssc, 1);
706 }
707
708 void CINTgout2e_int3c2e_spsp1(double *g,
709 double *gout, FINT *idx, CINTEnvVars *envs, FINT gout_empty);
int3c2e_spsp1_spinor_ssc(double complex * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)710 CACHE_SIZE_T int3c2e_spsp1_spinor_ssc(double complex *out, FINT *dims, FINT *shls, FINT *atm, FINT natm,
711 FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache)
712 {
713 FINT ng[] = {1, 1, 0, 0, 2, 4, 1, 1};
714 CINTEnvVars envs;
715 CINTinit_int3c2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
716 envs.f_gout = &CINTgout2e_int3c2e_spsp1;
717 return CINT3c2e_spinor_drv(out, dims, &envs, opt, cache, &c2s_si_3c2e1_ssc, 1);
718 }
int3c2e_ssc_optimizer(CINTOpt ** opt,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env)719 void int3c2e_ssc_optimizer(CINTOpt **opt, FINT *atm, FINT natm,
720 FINT *bas, FINT nbas, double *env)
721 {
722 int3c2e_optimizer(opt, atm, natm, bas, nbas, env);
723 }
724
725
726
727 ALL_CINT(int3c2e)
728 ALL_CINT_FORTRAN_(int3c2e)
729
730