1 #include "simint/boys/boys.h"
2 #include "simint/ostei/gen/ostei_generated.h"
3 #include "simint/vectorization/vectorization.h"
4 #include <math.h>
5 #include <string.h>
6 
7 
ostei_i_f_s_s(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_f_s_s)8 int ostei_i_f_s_s(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__i_f_s_s)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__i_f_s_s);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int iket;
26 
27     // partition workspace
28     double * const INT__i_s_s_s = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__k_s_s_s = work + (SIMINT_NSHELL_SIMD * 28);
30     double * const INT__l_s_s_s = work + (SIMINT_NSHELL_SIMD * 64);
31     double * const INT__m_s_s_s = work + (SIMINT_NSHELL_SIMD * 109);
32     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*164);
33     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
34     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 10;
35     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 37;
36     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 85;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 155;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 245;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 350;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 462;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 570;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_s_s = primwork + 660;
43     double * const hrrwork = (double *)(primwork + 715);
44     double * const HRR_INT__i_p_s_s = hrrwork + 0;
45     double * const HRR_INT__i_d_s_s = hrrwork + 84;
46     double * const HRR_INT__k_p_s_s = hrrwork + 252;
47     double * const HRR_INT__k_d_s_s = hrrwork + 360;
48     double * const HRR_INT__l_p_s_s = hrrwork + 576;
49 
50 
51     // Create constants
52     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
53     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
54     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
55     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
56     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
57     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
58     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
59     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
60     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
61 
62 
63     ////////////////////////////////////////
64     // Loop over shells and primitives
65     ////////////////////////////////////////
66 
67     real_abcd = 0;
68     istart = 0;
69     for(ab = 0; ab < P.nshell12_clip; ++ab)
70     {
71         const int iend = istart + P.nprim12[ab];
72 
73         cd = 0;
74         jstart = 0;
75 
76         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
77         {
78             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
79             int jend = jstart;
80             for(i = 0; i < nshellbatch; i++)
81                 jend += Q.nprim12[cd+i];
82 
83             // Clear the beginning of the workspace (where we are accumulating integrals)
84             memset(work, 0, SIMINT_NSHELL_SIMD * 164 * sizeof(double));
85             abcd = 0;
86 
87 
88             for(i = istart; i < iend; ++i)
89             {
90                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
91 
92                 if(check_screen)
93                 {
94                     // Skip this whole thing if always insignificant
95                     if((P.screen[i] * Q.screen_max) < screen_tol)
96                         continue;
97                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
98                 }
99 
100                 icd = 0;
101                 iprimcd = 0;
102                 nprim_icd = Q.nprim12[cd];
103                 double * restrict PRIM_PTR_INT__i_s_s_s = INT__i_s_s_s + abcd * 28;
104                 double * restrict PRIM_PTR_INT__k_s_s_s = INT__k_s_s_s + abcd * 36;
105                 double * restrict PRIM_PTR_INT__l_s_s_s = INT__l_s_s_s + abcd * 45;
106                 double * restrict PRIM_PTR_INT__m_s_s_s = INT__m_s_s_s + abcd * 55;
107 
108 
109 
110                 // Load these one per loop over i
111                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
112                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
113                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
114 
115                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
116 
117                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
118                 {
119                     // calculate the shell offsets
120                     // these are the offset from the shell pointed to by cd
121                     // for each element
122                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
123                     int lastoffset = 0;
124                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
125 
126                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
127                     {
128                         // Handle if the first element of the vector is a new shell
129                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
130                         {
131                             nprim_icd += Q.nprim12[cd + (++icd)];
132                             PRIM_PTR_INT__i_s_s_s += 28;
133                             PRIM_PTR_INT__k_s_s_s += 36;
134                             PRIM_PTR_INT__l_s_s_s += 45;
135                             PRIM_PTR_INT__m_s_s_s += 55;
136                         }
137                         iprimcd++;
138                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
139                         {
140                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
141                             {
142                                 shelloffsets[n] = shelloffsets[n-1] + 1;
143                                 lastoffset++;
144                                 nprim_icd += Q.nprim12[cd + (++icd)];
145                             }
146                             else
147                                 shelloffsets[n] = shelloffsets[n-1];
148                             iprimcd++;
149                         }
150                     }
151                     else
152                         iprimcd += SIMINT_SIMD_LEN;
153 
154                     // Do we have to compute this vector (or has it been screened out)?
155                     // (not_screened != 0 means we have to do this vector)
156                     if(check_screen)
157                     {
158                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
159                         if(vmax < screen_tol)
160                         {
161                             PRIM_PTR_INT__i_s_s_s += lastoffset*28;
162                             PRIM_PTR_INT__k_s_s_s += lastoffset*36;
163                             PRIM_PTR_INT__l_s_s_s += lastoffset*45;
164                             PRIM_PTR_INT__m_s_s_s += lastoffset*55;
165                             continue;
166                         }
167                     }
168 
169                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
170                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
171                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
172                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
173 
174 
175                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
176                     SIMINT_DBLTYPE PQ[3];
177                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
178                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
179                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
180                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
181                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
182                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
183 
184                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
185                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
186                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
187                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
188                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
189                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
190 
191                     // NOTE: Minus sign!
192                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
193                     SIMINT_DBLTYPE aop_PQ[3];
194                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
195                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
196                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
197 
198 
199                     //////////////////////////////////////////////
200                     // Fjt function section
201                     // Maximum v value: 9
202                     //////////////////////////////////////////////
203                     // The parameter to the Fjt function
204                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
205 
206 
207                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
208 
209 
210                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 9);
211                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
212                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
213                     for(n = 0; n <= 9; n++)
214                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
215 
216                     //////////////////////////////////////////////
217                     // Primitive integrals: Vertical recurrance
218                     //////////////////////////////////////////////
219 
220                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
221                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
222                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
223                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
224                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
225                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
226                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
227                     const SIMINT_DBLTYPE vrr_const_8_over_2p = SIMINT_MUL(const_8, one_over_2p);
228 
229 
230 
231                     // Forming PRIM_INT__p_s_s_s[9 * 3];
232                     for(n = 0; n < 9; ++n)  // loop over orders of auxiliary function
233                     {
234 
235                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
236                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]);
237 
238                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
239                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 1]);
240 
241                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
242                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 2]);
243 
244                     }
245 
246 
247 
248                     // Forming PRIM_INT__d_s_s_s[8 * 6];
249                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
250                     {
251 
252                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
253                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 0]);
254                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 0]);
255 
256                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
257                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 3]);
258                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 3]);
259 
260                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
261                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_s_s[n * 6 + 5]);
262                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 5]);
263 
264                     }
265 
266 
267 
268                     // Forming PRIM_INT__f_s_s_s[7 * 10];
269                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
270                     {
271 
272                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
273                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 0]);
274                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__f_s_s_s[n * 10 + 0]);
275 
276                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
277                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 1]);
278 
279                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
280                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 6]);
281                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__f_s_s_s[n * 10 + 6]);
282 
283                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
284                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 9]);
285                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__f_s_s_s[n * 10 + 9]);
286 
287                     }
288 
289 
290                     VRR_I_g_s_s_s(
291                             PRIM_INT__g_s_s_s,
292                             PRIM_INT__f_s_s_s,
293                             PRIM_INT__d_s_s_s,
294                             P_PA,
295                             a_over_p,
296                             aop_PQ,
297                             one_over_2p,
298                             6);
299 
300 
301                     VRR_I_h_s_s_s(
302                             PRIM_INT__h_s_s_s,
303                             PRIM_INT__g_s_s_s,
304                             PRIM_INT__f_s_s_s,
305                             P_PA,
306                             a_over_p,
307                             aop_PQ,
308                             one_over_2p,
309                             5);
310 
311 
312                     ostei_general_vrr1_I(6, 4,
313                             one_over_2p, a_over_p, aop_PQ, P_PA,
314                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
315 
316 
317                     ostei_general_vrr1_I(7, 3,
318                             one_over_2p, a_over_p, aop_PQ, P_PA,
319                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
320 
321 
322                     ostei_general_vrr1_I(8, 2,
323                             one_over_2p, a_over_p, aop_PQ, P_PA,
324                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
325 
326 
327                     ostei_general_vrr1_I(9, 1,
328                             one_over_2p, a_over_p, aop_PQ, P_PA,
329                             PRIM_INT__l_s_s_s, PRIM_INT__k_s_s_s, PRIM_INT__m_s_s_s);
330 
331 
332 
333 
334                     ////////////////////////////////////
335                     // Accumulate contracted integrals
336                     ////////////////////////////////////
337                     if(lastoffset == 0)
338                     {
339                         contract_all(28, PRIM_INT__i_s_s_s, PRIM_PTR_INT__i_s_s_s);
340                         contract_all(36, PRIM_INT__k_s_s_s, PRIM_PTR_INT__k_s_s_s);
341                         contract_all(45, PRIM_INT__l_s_s_s, PRIM_PTR_INT__l_s_s_s);
342                         contract_all(55, PRIM_INT__m_s_s_s, PRIM_PTR_INT__m_s_s_s);
343                     }
344                     else
345                     {
346                         contract(28, shelloffsets, PRIM_INT__i_s_s_s, PRIM_PTR_INT__i_s_s_s);
347                         contract(36, shelloffsets, PRIM_INT__k_s_s_s, PRIM_PTR_INT__k_s_s_s);
348                         contract(45, shelloffsets, PRIM_INT__l_s_s_s, PRIM_PTR_INT__l_s_s_s);
349                         contract(55, shelloffsets, PRIM_INT__m_s_s_s, PRIM_PTR_INT__m_s_s_s);
350                         PRIM_PTR_INT__i_s_s_s += lastoffset*28;
351                         PRIM_PTR_INT__k_s_s_s += lastoffset*36;
352                         PRIM_PTR_INT__l_s_s_s += lastoffset*45;
353                         PRIM_PTR_INT__m_s_s_s += lastoffset*55;
354                     }
355 
356                 }  // close loop over j
357             }  // close loop over i
358 
359             //Advance to the next batch
360             jstart = SIMINT_SIMD_ROUND(jend);
361 
362             //////////////////////////////////////////////
363             // Contracted integrals: Horizontal recurrance
364             //////////////////////////////////////////////
365 
366 
367             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
368 
369 
370             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
371             {
372 
373                 // set up HRR pointers
374                 double const * restrict HRR_INT__i_s_s_s = INT__i_s_s_s + abcd * 28;
375                 double const * restrict HRR_INT__k_s_s_s = INT__k_s_s_s + abcd * 36;
376                 double const * restrict HRR_INT__l_s_s_s = INT__l_s_s_s + abcd * 45;
377                 double const * restrict HRR_INT__m_s_s_s = INT__m_s_s_s + abcd * 55;
378                 double * restrict HRR_INT__i_f_s_s = INT__i_f_s_s + real_abcd * 280;
379 
380                 // form INT__i_p_s_s
381                 ostei_general_hrr_J(6, 1, 0, 0, hAB, HRR_INT__k_s_s_s, HRR_INT__i_s_s_s, HRR_INT__i_p_s_s);
382 
383                 // form INT__k_p_s_s
384                 ostei_general_hrr_J(7, 1, 0, 0, hAB, HRR_INT__l_s_s_s, HRR_INT__k_s_s_s, HRR_INT__k_p_s_s);
385 
386                 // form INT__l_p_s_s
387                 ostei_general_hrr_J(8, 1, 0, 0, hAB, HRR_INT__m_s_s_s, HRR_INT__l_s_s_s, HRR_INT__l_p_s_s);
388 
389                 // form INT__i_d_s_s
390                 ostei_general_hrr_J(6, 2, 0, 0, hAB, HRR_INT__k_p_s_s, HRR_INT__i_p_s_s, HRR_INT__i_d_s_s);
391 
392                 // form INT__k_d_s_s
393                 ostei_general_hrr_J(7, 2, 0, 0, hAB, HRR_INT__l_p_s_s, HRR_INT__k_p_s_s, HRR_INT__k_d_s_s);
394 
395                 // form INT__i_f_s_s
396                 ostei_general_hrr_J(6, 3, 0, 0, hAB, HRR_INT__k_d_s_s, HRR_INT__i_d_s_s, HRR_INT__i_f_s_s);
397 
398 
399             }  // close HRR loop
400 
401 
402         }   // close loop cdbatch
403 
404         istart = iend;
405     }  // close loop over ab
406 
407     return P.nshell12_clip * Q.nshell12_clip;
408 }
409 
ostei_f_i_s_s(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__f_i_s_s)410 int ostei_f_i_s_s(struct simint_multi_shellpair const P,
411                   struct simint_multi_shellpair const Q,
412                   double screen_tol,
413                   double * const restrict work,
414                   double * const restrict INT__f_i_s_s)
415 {
416     double P_AB[3*P.nshell12];
417     struct simint_multi_shellpair P_tmp = P;
418     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
419     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
420     P_tmp.AB_x = P_AB;
421     P_tmp.AB_y = P_AB + P.nshell12;
422     P_tmp.AB_z = P_AB + 2*P.nshell12;
423 
424     for(int i = 0; i < P.nshell12; i++)
425     {
426         P_tmp.AB_x[i] = -P.AB_x[i];
427         P_tmp.AB_y[i] = -P.AB_y[i];
428         P_tmp.AB_z[i] = -P.AB_z[i];
429     }
430 
431     int ret = ostei_i_f_s_s(P_tmp, Q, screen_tol, work, INT__f_i_s_s);
432     double buffer[280] SIMINT_ALIGN_ARRAY_DBL;
433 
434     for(int q = 0; q < ret; q++)
435     {
436         int idx = 0;
437         for(int a = 0; a < 10; ++a)
438         for(int b = 0; b < 28; ++b)
439         for(int c = 0; c < 1; ++c)
440         for(int d = 0; d < 1; ++d)
441             buffer[idx++] = INT__f_i_s_s[q*280+b*10+a*1+c*1+d];
442 
443         memcpy(INT__f_i_s_s+q*280, buffer, 280*sizeof(double));
444     }
445 
446     return ret;
447 }
448 
449