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_p_f_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_p_f_d)8 int ostei_i_p_f_d(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_p_f_d)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__i_p_f_d);
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     int ibra;
27 
28     // partition workspace
29     double * const INT__i_s_f_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__i_s_g_s = work + (SIMINT_NSHELL_SIMD * 280);
31     double * const INT__i_s_h_s = work + (SIMINT_NSHELL_SIMD * 700);
32     double * const INT__k_s_f_s = work + (SIMINT_NSHELL_SIMD * 1288);
33     double * const INT__k_s_g_s = work + (SIMINT_NSHELL_SIMD * 1648);
34     double * const INT__k_s_h_s = work + (SIMINT_NSHELL_SIMD * 2188);
35     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*2944);
36     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 13;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 49;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 115;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 205;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 305;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 455;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 695;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 830;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 1055;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 1415;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 1865;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 2033;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 2348;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 2852;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_g_s = primwork + 3482;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 4112;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 4308;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 4728;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 5400;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 6240;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_h_s = primwork + 7080;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 7668;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 7884;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 8424;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 9288;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_g_s = primwork + 10368;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_h_s = primwork + 11448;
64     double * const hrrwork = (double *)(primwork + 12204);
65     double * const HRR_INT__i_p_f_s = hrrwork + 0;
66     double * const HRR_INT__i_p_f_p = hrrwork + 840;
67     double * const HRR_INT__i_p_g_s = hrrwork + 3360;
68     double * const HRR_INT__i_p_g_p = hrrwork + 4620;
69     double * const HRR_INT__i_p_h_s = hrrwork + 8400;
70 
71 
72     // Create constants
73     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
74     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
75     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
76     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
77     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
78     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
79     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
80     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
81 
82 
83     ////////////////////////////////////////
84     // Loop over shells and primitives
85     ////////////////////////////////////////
86 
87     real_abcd = 0;
88     istart = 0;
89     for(ab = 0; ab < P.nshell12_clip; ++ab)
90     {
91         const int iend = istart + P.nprim12[ab];
92 
93         cd = 0;
94         jstart = 0;
95 
96         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
97         {
98             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
99             int jend = jstart;
100             for(i = 0; i < nshellbatch; i++)
101                 jend += Q.nprim12[cd+i];
102 
103             // Clear the beginning of the workspace (where we are accumulating integrals)
104             memset(work, 0, SIMINT_NSHELL_SIMD * 2944 * sizeof(double));
105             abcd = 0;
106 
107 
108             for(i = istart; i < iend; ++i)
109             {
110                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
111 
112                 if(check_screen)
113                 {
114                     // Skip this whole thing if always insignificant
115                     if((P.screen[i] * Q.screen_max) < screen_tol)
116                         continue;
117                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
118                 }
119 
120                 icd = 0;
121                 iprimcd = 0;
122                 nprim_icd = Q.nprim12[cd];
123                 double * restrict PRIM_PTR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
124                 double * restrict PRIM_PTR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
125                 double * restrict PRIM_PTR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
126                 double * restrict PRIM_PTR_INT__k_s_f_s = INT__k_s_f_s + abcd * 360;
127                 double * restrict PRIM_PTR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
128                 double * restrict PRIM_PTR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
129 
130 
131 
132                 // Load these one per loop over i
133                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
134                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
135                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
136 
137                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
138 
139                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
140                 {
141                     // calculate the shell offsets
142                     // these are the offset from the shell pointed to by cd
143                     // for each element
144                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
145                     int lastoffset = 0;
146                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
147 
148                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
149                     {
150                         // Handle if the first element of the vector is a new shell
151                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
152                         {
153                             nprim_icd += Q.nprim12[cd + (++icd)];
154                             PRIM_PTR_INT__i_s_f_s += 280;
155                             PRIM_PTR_INT__i_s_g_s += 420;
156                             PRIM_PTR_INT__i_s_h_s += 588;
157                             PRIM_PTR_INT__k_s_f_s += 360;
158                             PRIM_PTR_INT__k_s_g_s += 540;
159                             PRIM_PTR_INT__k_s_h_s += 756;
160                         }
161                         iprimcd++;
162                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
163                         {
164                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
165                             {
166                                 shelloffsets[n] = shelloffsets[n-1] + 1;
167                                 lastoffset++;
168                                 nprim_icd += Q.nprim12[cd + (++icd)];
169                             }
170                             else
171                                 shelloffsets[n] = shelloffsets[n-1];
172                             iprimcd++;
173                         }
174                     }
175                     else
176                         iprimcd += SIMINT_SIMD_LEN;
177 
178                     // Do we have to compute this vector (or has it been screened out)?
179                     // (not_screened != 0 means we have to do this vector)
180                     if(check_screen)
181                     {
182                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
183                         if(vmax < screen_tol)
184                         {
185                             PRIM_PTR_INT__i_s_f_s += lastoffset*280;
186                             PRIM_PTR_INT__i_s_g_s += lastoffset*420;
187                             PRIM_PTR_INT__i_s_h_s += lastoffset*588;
188                             PRIM_PTR_INT__k_s_f_s += lastoffset*360;
189                             PRIM_PTR_INT__k_s_g_s += lastoffset*540;
190                             PRIM_PTR_INT__k_s_h_s += lastoffset*756;
191                             continue;
192                         }
193                     }
194 
195                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
196                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
197                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
198                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
199 
200 
201                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
202                     SIMINT_DBLTYPE PQ[3];
203                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
204                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
205                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
206                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
207                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
208                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
209 
210                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
211                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
212                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
213                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
214                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
215                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
216                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
217 
218                     // NOTE: Minus sign!
219                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
220                     SIMINT_DBLTYPE aop_PQ[3];
221                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
222                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
223                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
224 
225                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
226                     SIMINT_DBLTYPE aoq_PQ[3];
227                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
228                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
229                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
230                     // Put a minus sign here so we don't have to in RR routines
231                     a_over_q = SIMINT_NEG(a_over_q);
232 
233 
234                     //////////////////////////////////////////////
235                     // Fjt function section
236                     // Maximum v value: 12
237                     //////////////////////////////////////////////
238                     // The parameter to the Fjt function
239                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
240 
241 
242                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
243 
244 
245                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 12);
246                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
247                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
248                     for(n = 0; n <= 12; n++)
249                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
250 
251                     //////////////////////////////////////////////
252                     // Primitive integrals: Vertical recurrance
253                     //////////////////////////////////////////////
254 
255                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
256                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
257                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
258                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
259                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
260                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
261                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
262                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
263                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
264                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
265                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
266                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
267                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
268                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
269                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
270                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
271                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
272 
273 
274 
275                     // Forming PRIM_INT__p_s_s_s[12 * 3];
276                     for(n = 0; n < 12; ++n)  // loop over orders of auxiliary function
277                     {
278 
279                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
280                         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]);
281 
282                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
283                         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]);
284 
285                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
286                         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]);
287 
288                     }
289 
290 
291 
292                     // Forming PRIM_INT__d_s_s_s[11 * 6];
293                     for(n = 0; n < 11; ++n)  // loop over orders of auxiliary function
294                     {
295 
296                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
297                         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]);
298                         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]);
299 
300                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
301                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 1]);
302 
303                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
304                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 2]);
305 
306                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
307                         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]);
308                         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]);
309 
310                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
311                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 4]);
312 
313                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
314                         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]);
315                         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]);
316 
317                     }
318 
319 
320 
321                     // Forming PRIM_INT__f_s_s_s[10 * 10];
322                     for(n = 0; n < 10; ++n)  // loop over orders of auxiliary function
323                     {
324 
325                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
326                         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]);
327                         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]);
328 
329                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
330                         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]);
331 
332                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
333                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 2]);
334 
335                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
336                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 3]);
337 
338                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
339                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__f_s_s_s[n * 10 + 4]);
340 
341                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
342                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 5]);
343 
344                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
345                         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]);
346                         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]);
347 
348                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
349                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 7]);
350 
351                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
352                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 8]);
353 
354                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
355                         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]);
356                         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]);
357 
358                     }
359 
360 
361                     VRR_I_g_s_s_s(
362                             PRIM_INT__g_s_s_s,
363                             PRIM_INT__f_s_s_s,
364                             PRIM_INT__d_s_s_s,
365                             P_PA,
366                             a_over_p,
367                             aop_PQ,
368                             one_over_2p,
369                             9);
370 
371 
372                     VRR_I_h_s_s_s(
373                             PRIM_INT__h_s_s_s,
374                             PRIM_INT__g_s_s_s,
375                             PRIM_INT__f_s_s_s,
376                             P_PA,
377                             a_over_p,
378                             aop_PQ,
379                             one_over_2p,
380                             8);
381 
382 
383                     ostei_general_vrr1_I(6, 7,
384                             one_over_2p, a_over_p, aop_PQ, P_PA,
385                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
386 
387 
388                     ostei_general_vrr_K(6, 0, 1, 0, 5,
389                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
390                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
391 
392 
393                     ostei_general_vrr_K(5, 0, 1, 0, 5,
394                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
395                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
396 
397 
398                     ostei_general_vrr_K(6, 0, 2, 0, 4,
399                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
400                             PRIM_INT__i_s_p_s, PRIM_INT__i_s_s_s, NULL, PRIM_INT__h_s_p_s, NULL, PRIM_INT__i_s_d_s);
401 
402 
403                     VRR_K_g_s_p_s(
404                             PRIM_INT__g_s_p_s,
405                             PRIM_INT__g_s_s_s,
406                             PRIM_INT__f_s_s_s,
407                             Q_PA,
408                             aoq_PQ,
409                             one_over_2pq,
410                             5);
411 
412 
413                     ostei_general_vrr_K(5, 0, 2, 0, 4,
414                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
415                             PRIM_INT__h_s_p_s, PRIM_INT__h_s_s_s, NULL, PRIM_INT__g_s_p_s, NULL, PRIM_INT__h_s_d_s);
416 
417 
418                     ostei_general_vrr_K(6, 0, 3, 0, 3,
419                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
420                             PRIM_INT__i_s_d_s, PRIM_INT__i_s_p_s, NULL, PRIM_INT__h_s_d_s, NULL, PRIM_INT__i_s_f_s);
421 
422 
423                     VRR_K_f_s_p_s(
424                             PRIM_INT__f_s_p_s,
425                             PRIM_INT__f_s_s_s,
426                             PRIM_INT__d_s_s_s,
427                             Q_PA,
428                             aoq_PQ,
429                             one_over_2pq,
430                             5);
431 
432 
433                     ostei_general_vrr_K(4, 0, 2, 0, 4,
434                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
435                             PRIM_INT__g_s_p_s, PRIM_INT__g_s_s_s, NULL, PRIM_INT__f_s_p_s, NULL, PRIM_INT__g_s_d_s);
436 
437 
438                     ostei_general_vrr_K(5, 0, 3, 0, 3,
439                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
440                             PRIM_INT__h_s_d_s, PRIM_INT__h_s_p_s, NULL, PRIM_INT__g_s_d_s, NULL, PRIM_INT__h_s_f_s);
441 
442 
443                     ostei_general_vrr_K(6, 0, 4, 0, 2,
444                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
445                             PRIM_INT__i_s_f_s, PRIM_INT__i_s_d_s, NULL, PRIM_INT__h_s_f_s, NULL, PRIM_INT__i_s_g_s);
446 
447 
448                     ostei_general_vrr1_I(7, 6,
449                             one_over_2p, a_over_p, aop_PQ, P_PA,
450                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
451 
452 
453                     ostei_general_vrr_K(7, 0, 1, 0, 5,
454                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
455                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
456 
457 
458                     ostei_general_vrr_K(7, 0, 2, 0, 4,
459                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
460                             PRIM_INT__k_s_p_s, PRIM_INT__k_s_s_s, NULL, PRIM_INT__i_s_p_s, NULL, PRIM_INT__k_s_d_s);
461 
462 
463                     ostei_general_vrr_K(7, 0, 3, 0, 3,
464                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
465                             PRIM_INT__k_s_d_s, PRIM_INT__k_s_p_s, NULL, PRIM_INT__i_s_d_s, NULL, PRIM_INT__k_s_f_s);
466 
467 
468 
469                     // Forming PRIM_INT__d_s_p_s[5 * 18];
470                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
471                     {
472 
473                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
474                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
475                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
476 
477                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
478                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 1]);
479 
480                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
481                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 2]);
482 
483                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
484                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
485                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
486 
487                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
488                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 4]);
489                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 4]);
490 
491                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
492                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 5]);
493 
494                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
495                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
496                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
497 
498                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
499                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 7]);
500 
501                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
502                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 8]);
503                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 8]);
504 
505                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
506                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 9]);
507 
508                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
509                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 10]);
510                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 10]);
511 
512                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
513                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 11]);
514 
515                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
516                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 12]);
517 
518                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
519                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 13]);
520                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 13]);
521 
522                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
523                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 14]);
524                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 14]);
525 
526                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
527                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 15]);
528 
529                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
530                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 16]);
531 
532                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
533                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 17]);
534                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 17]);
535 
536                     }
537 
538 
539                     VRR_K_f_s_d_s(
540                             PRIM_INT__f_s_d_s,
541                             PRIM_INT__f_s_p_s,
542                             PRIM_INT__f_s_s_s,
543                             PRIM_INT__d_s_p_s,
544                             Q_PA,
545                             a_over_q,
546                             aoq_PQ,
547                             one_over_2pq,
548                             one_over_2q,
549                             4);
550 
551 
552                     ostei_general_vrr_K(4, 0, 3, 0, 3,
553                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
554                             PRIM_INT__g_s_d_s, PRIM_INT__g_s_p_s, NULL, PRIM_INT__f_s_d_s, NULL, PRIM_INT__g_s_f_s);
555 
556 
557                     ostei_general_vrr_K(5, 0, 4, 0, 2,
558                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
559                             PRIM_INT__h_s_f_s, PRIM_INT__h_s_d_s, NULL, PRIM_INT__g_s_f_s, NULL, PRIM_INT__h_s_g_s);
560 
561 
562                     ostei_general_vrr_K(6, 0, 5, 0, 1,
563                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
564                             PRIM_INT__i_s_g_s, PRIM_INT__i_s_f_s, NULL, PRIM_INT__h_s_g_s, NULL, PRIM_INT__i_s_h_s);
565 
566 
567                     ostei_general_vrr_K(7, 0, 4, 0, 2,
568                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
569                             PRIM_INT__k_s_f_s, PRIM_INT__k_s_d_s, NULL, PRIM_INT__i_s_f_s, NULL, PRIM_INT__k_s_g_s);
570 
571 
572                     ostei_general_vrr_K(7, 0, 5, 0, 1,
573                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
574                             PRIM_INT__k_s_g_s, PRIM_INT__k_s_f_s, NULL, PRIM_INT__i_s_g_s, NULL, PRIM_INT__k_s_h_s);
575 
576 
577 
578 
579                     ////////////////////////////////////
580                     // Accumulate contracted integrals
581                     ////////////////////////////////////
582                     if(lastoffset == 0)
583                     {
584                         contract_all(280, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
585                         contract_all(420, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
586                         contract_all(588, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
587                         contract_all(360, PRIM_INT__k_s_f_s, PRIM_PTR_INT__k_s_f_s);
588                         contract_all(540, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
589                         contract_all(756, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
590                     }
591                     else
592                     {
593                         contract(280, shelloffsets, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
594                         contract(420, shelloffsets, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
595                         contract(588, shelloffsets, PRIM_INT__i_s_h_s, PRIM_PTR_INT__i_s_h_s);
596                         contract(360, shelloffsets, PRIM_INT__k_s_f_s, PRIM_PTR_INT__k_s_f_s);
597                         contract(540, shelloffsets, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
598                         contract(756, shelloffsets, PRIM_INT__k_s_h_s, PRIM_PTR_INT__k_s_h_s);
599                         PRIM_PTR_INT__i_s_f_s += lastoffset*280;
600                         PRIM_PTR_INT__i_s_g_s += lastoffset*420;
601                         PRIM_PTR_INT__i_s_h_s += lastoffset*588;
602                         PRIM_PTR_INT__k_s_f_s += lastoffset*360;
603                         PRIM_PTR_INT__k_s_g_s += lastoffset*540;
604                         PRIM_PTR_INT__k_s_h_s += lastoffset*756;
605                     }
606 
607                 }  // close loop over j
608             }  // close loop over i
609 
610             //Advance to the next batch
611             jstart = SIMINT_SIMD_ROUND(jend);
612 
613             //////////////////////////////////////////////
614             // Contracted integrals: Horizontal recurrance
615             //////////////////////////////////////////////
616 
617 
618             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
619 
620 
621             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
622             {
623                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
624 
625                 // set up HRR pointers
626                 double const * restrict HRR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
627                 double const * restrict HRR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
628                 double const * restrict HRR_INT__i_s_h_s = INT__i_s_h_s + abcd * 588;
629                 double const * restrict HRR_INT__k_s_f_s = INT__k_s_f_s + abcd * 360;
630                 double const * restrict HRR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
631                 double const * restrict HRR_INT__k_s_h_s = INT__k_s_h_s + abcd * 756;
632                 double * restrict HRR_INT__i_p_f_d = INT__i_p_f_d + real_abcd * 5040;
633 
634                 // form INT__i_p_f_s
635                 ostei_general_hrr_J(6, 1, 3, 0, hAB, HRR_INT__k_s_f_s, HRR_INT__i_s_f_s, HRR_INT__i_p_f_s);
636 
637                 // form INT__i_p_g_s
638                 ostei_general_hrr_J(6, 1, 4, 0, hAB, HRR_INT__k_s_g_s, HRR_INT__i_s_g_s, HRR_INT__i_p_g_s);
639 
640                 // form INT__i_p_h_s
641                 ostei_general_hrr_J(6, 1, 5, 0, hAB, HRR_INT__k_s_h_s, HRR_INT__i_s_h_s, HRR_INT__i_p_h_s);
642 
643                 // form INT__i_p_f_p
644                 HRR_L_f_p(
645                     HRR_INT__i_p_f_p,
646                     HRR_INT__i_p_f_s,
647                     HRR_INT__i_p_g_s,
648                     hCD, 84);
649 
650                 // form INT__i_p_g_p
651                 HRR_L_g_p(
652                     HRR_INT__i_p_g_p,
653                     HRR_INT__i_p_g_s,
654                     HRR_INT__i_p_h_s,
655                     hCD, 84);
656 
657                 // form INT__i_p_f_d
658                 HRR_L_f_d(
659                     HRR_INT__i_p_f_d,
660                     HRR_INT__i_p_f_p,
661                     HRR_INT__i_p_g_p,
662                     hCD, 84);
663 
664 
665             }  // close HRR loop
666 
667 
668         }   // close loop cdbatch
669 
670         istart = iend;
671     }  // close loop over ab
672 
673     return P.nshell12_clip * Q.nshell12_clip;
674 }
675 
ostei_p_i_f_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__p_i_f_d)676 int ostei_p_i_f_d(struct simint_multi_shellpair const P,
677                   struct simint_multi_shellpair const Q,
678                   double screen_tol,
679                   double * const restrict work,
680                   double * const restrict INT__p_i_f_d)
681 {
682     double P_AB[3*P.nshell12];
683     struct simint_multi_shellpair P_tmp = P;
684     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
685     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
686     P_tmp.AB_x = P_AB;
687     P_tmp.AB_y = P_AB + P.nshell12;
688     P_tmp.AB_z = P_AB + 2*P.nshell12;
689 
690     for(int i = 0; i < P.nshell12; i++)
691     {
692         P_tmp.AB_x[i] = -P.AB_x[i];
693         P_tmp.AB_y[i] = -P.AB_y[i];
694         P_tmp.AB_z[i] = -P.AB_z[i];
695     }
696 
697     int ret = ostei_i_p_f_d(P_tmp, Q, screen_tol, work, INT__p_i_f_d);
698     double buffer[5040] SIMINT_ALIGN_ARRAY_DBL;
699 
700     for(int q = 0; q < ret; q++)
701     {
702         int idx = 0;
703         for(int a = 0; a < 3; ++a)
704         for(int b = 0; b < 28; ++b)
705         for(int c = 0; c < 10; ++c)
706         for(int d = 0; d < 6; ++d)
707             buffer[idx++] = INT__p_i_f_d[q*5040+b*180+a*60+c*6+d];
708 
709         memcpy(INT__p_i_f_d+q*5040, buffer, 5040*sizeof(double));
710     }
711 
712     return ret;
713 }
714 
ostei_i_p_d_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_p_d_f)715 int ostei_i_p_d_f(struct simint_multi_shellpair const P,
716                   struct simint_multi_shellpair const Q,
717                   double screen_tol,
718                   double * const restrict work,
719                   double * const restrict INT__i_p_d_f)
720 {
721     double Q_AB[3*Q.nshell12];
722     struct simint_multi_shellpair Q_tmp = Q;
723     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
724     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
725     Q_tmp.AB_x = Q_AB;
726     Q_tmp.AB_y = Q_AB + Q.nshell12;
727     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
728 
729     for(int i = 0; i < Q.nshell12; i++)
730     {
731         Q_tmp.AB_x[i] = -Q.AB_x[i];
732         Q_tmp.AB_y[i] = -Q.AB_y[i];
733         Q_tmp.AB_z[i] = -Q.AB_z[i];
734     }
735 
736     int ret = ostei_i_p_f_d(P, Q_tmp, screen_tol, work, INT__i_p_d_f);
737     double buffer[5040] SIMINT_ALIGN_ARRAY_DBL;
738 
739     for(int q = 0; q < ret; q++)
740     {
741         int idx = 0;
742         for(int a = 0; a < 28; ++a)
743         for(int b = 0; b < 3; ++b)
744         for(int c = 0; c < 6; ++c)
745         for(int d = 0; d < 10; ++d)
746             buffer[idx++] = INT__i_p_d_f[q*5040+a*180+b*60+d*6+c];
747 
748         memcpy(INT__i_p_d_f+q*5040, buffer, 5040*sizeof(double));
749     }
750 
751     return ret;
752 }
753 
ostei_p_i_d_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__p_i_d_f)754 int ostei_p_i_d_f(struct simint_multi_shellpair const P,
755                   struct simint_multi_shellpair const Q,
756                   double screen_tol,
757                   double * const restrict work,
758                   double * const restrict INT__p_i_d_f)
759 {
760     double P_AB[3*P.nshell12];
761     struct simint_multi_shellpair P_tmp = P;
762     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
763     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
764     P_tmp.AB_x = P_AB;
765     P_tmp.AB_y = P_AB + P.nshell12;
766     P_tmp.AB_z = P_AB + 2*P.nshell12;
767 
768     for(int i = 0; i < P.nshell12; i++)
769     {
770         P_tmp.AB_x[i] = -P.AB_x[i];
771         P_tmp.AB_y[i] = -P.AB_y[i];
772         P_tmp.AB_z[i] = -P.AB_z[i];
773     }
774 
775     double Q_AB[3*Q.nshell12];
776     struct simint_multi_shellpair Q_tmp = Q;
777     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
778     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
779     Q_tmp.AB_x = Q_AB;
780     Q_tmp.AB_y = Q_AB + Q.nshell12;
781     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
782 
783     for(int i = 0; i < Q.nshell12; i++)
784     {
785         Q_tmp.AB_x[i] = -Q.AB_x[i];
786         Q_tmp.AB_y[i] = -Q.AB_y[i];
787         Q_tmp.AB_z[i] = -Q.AB_z[i];
788     }
789 
790     int ret = ostei_i_p_f_d(P_tmp, Q_tmp, screen_tol, work, INT__p_i_d_f);
791     double buffer[5040] SIMINT_ALIGN_ARRAY_DBL;
792 
793     for(int q = 0; q < ret; q++)
794     {
795         int idx = 0;
796         for(int a = 0; a < 3; ++a)
797         for(int b = 0; b < 28; ++b)
798         for(int c = 0; c < 6; ++c)
799         for(int d = 0; d < 10; ++d)
800             buffer[idx++] = INT__p_i_d_f[q*5040+b*180+a*60+d*6+c];
801 
802         memcpy(INT__p_i_d_f+q*5040, buffer, 5040*sizeof(double));
803     }
804 
805     return ret;
806 }
807 
808