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_f_f_f_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_f_f_s)8 int ostei_f_f_f_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__f_f_f_s)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__f_f_f_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__f_s_f_s = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__g_s_f_s = work + (SIMINT_NSHELL_SIMD * 100);
30     double * const INT__h_s_f_s = work + (SIMINT_NSHELL_SIMD * 250);
31     double * const INT__i_s_f_s = work + (SIMINT_NSHELL_SIMD * 460);
32     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*740);
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__p_s_p_s = primwork + 37;
36     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 64;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 112;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_d_s = primwork + 166;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 238;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 308;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 398;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_f_s = primwork + 518;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 618;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 708;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 843;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 1023;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 1173;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 1278;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 1467;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 1719;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 1929;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 2041;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 2293;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 2629;
55     double * const hrrwork = (double *)(primwork + 2909);
56     double * const HRR_INT__f_p_f_s = hrrwork + 0;
57     double * const HRR_INT__f_d_f_s = hrrwork + 300;
58     double * const HRR_INT__g_p_f_s = hrrwork + 900;
59     double * const HRR_INT__g_d_f_s = hrrwork + 1350;
60     double * const HRR_INT__h_p_f_s = hrrwork + 2250;
61 
62 
63     // Create constants
64     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
65     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
66     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
67     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
68     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
69     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
70     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
71 
72 
73     ////////////////////////////////////////
74     // Loop over shells and primitives
75     ////////////////////////////////////////
76 
77     real_abcd = 0;
78     istart = 0;
79     for(ab = 0; ab < P.nshell12_clip; ++ab)
80     {
81         const int iend = istart + P.nprim12[ab];
82 
83         cd = 0;
84         jstart = 0;
85 
86         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
87         {
88             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
89             int jend = jstart;
90             for(i = 0; i < nshellbatch; i++)
91                 jend += Q.nprim12[cd+i];
92 
93             // Clear the beginning of the workspace (where we are accumulating integrals)
94             memset(work, 0, SIMINT_NSHELL_SIMD * 740 * sizeof(double));
95             abcd = 0;
96 
97 
98             for(i = istart; i < iend; ++i)
99             {
100                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
101 
102                 if(check_screen)
103                 {
104                     // Skip this whole thing if always insignificant
105                     if((P.screen[i] * Q.screen_max) < screen_tol)
106                         continue;
107                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
108                 }
109 
110                 icd = 0;
111                 iprimcd = 0;
112                 nprim_icd = Q.nprim12[cd];
113                 double * restrict PRIM_PTR_INT__f_s_f_s = INT__f_s_f_s + abcd * 100;
114                 double * restrict PRIM_PTR_INT__g_s_f_s = INT__g_s_f_s + abcd * 150;
115                 double * restrict PRIM_PTR_INT__h_s_f_s = INT__h_s_f_s + abcd * 210;
116                 double * restrict PRIM_PTR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
117 
118 
119 
120                 // Load these one per loop over i
121                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
122                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
123                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
124 
125                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
126 
127                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
128                 {
129                     // calculate the shell offsets
130                     // these are the offset from the shell pointed to by cd
131                     // for each element
132                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
133                     int lastoffset = 0;
134                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
135 
136                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
137                     {
138                         // Handle if the first element of the vector is a new shell
139                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
140                         {
141                             nprim_icd += Q.nprim12[cd + (++icd)];
142                             PRIM_PTR_INT__f_s_f_s += 100;
143                             PRIM_PTR_INT__g_s_f_s += 150;
144                             PRIM_PTR_INT__h_s_f_s += 210;
145                             PRIM_PTR_INT__i_s_f_s += 280;
146                         }
147                         iprimcd++;
148                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
149                         {
150                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
151                             {
152                                 shelloffsets[n] = shelloffsets[n-1] + 1;
153                                 lastoffset++;
154                                 nprim_icd += Q.nprim12[cd + (++icd)];
155                             }
156                             else
157                                 shelloffsets[n] = shelloffsets[n-1];
158                             iprimcd++;
159                         }
160                     }
161                     else
162                         iprimcd += SIMINT_SIMD_LEN;
163 
164                     // Do we have to compute this vector (or has it been screened out)?
165                     // (not_screened != 0 means we have to do this vector)
166                     if(check_screen)
167                     {
168                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
169                         if(vmax < screen_tol)
170                         {
171                             PRIM_PTR_INT__f_s_f_s += lastoffset*100;
172                             PRIM_PTR_INT__g_s_f_s += lastoffset*150;
173                             PRIM_PTR_INT__h_s_f_s += lastoffset*210;
174                             PRIM_PTR_INT__i_s_f_s += lastoffset*280;
175                             continue;
176                         }
177                     }
178 
179                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
180                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
181                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
182                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
183 
184 
185                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
186                     SIMINT_DBLTYPE PQ[3];
187                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
188                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
189                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
190                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
191                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
192                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
193 
194                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
195                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
196                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
197                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
198                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
199                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
200                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
201 
202                     // NOTE: Minus sign!
203                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
204                     SIMINT_DBLTYPE aop_PQ[3];
205                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
206                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
207                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
208 
209                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
210                     SIMINT_DBLTYPE aoq_PQ[3];
211                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
212                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
213                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
214                     // Put a minus sign here so we don't have to in RR routines
215                     a_over_q = SIMINT_NEG(a_over_q);
216 
217 
218                     //////////////////////////////////////////////
219                     // Fjt function section
220                     // Maximum v value: 9
221                     //////////////////////////////////////////////
222                     // The parameter to the Fjt function
223                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
224 
225 
226                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
227 
228 
229                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 9);
230                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
231                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
232                     for(n = 0; n <= 9; n++)
233                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
234 
235                     //////////////////////////////////////////////
236                     // Primitive integrals: Vertical recurrance
237                     //////////////////////////////////////////////
238 
239                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
240                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
241                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
242                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
243                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
244                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
245                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
246                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
247                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
248                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
249                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
250                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
251                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
252 
253 
254 
255                     // Forming PRIM_INT__p_s_s_s[9 * 3];
256                     for(n = 0; n < 9; ++n)  // loop over orders of auxiliary function
257                     {
258 
259                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
260                         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]);
261 
262                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
263                         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]);
264 
265                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
266                         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]);
267 
268                     }
269 
270 
271 
272                     // Forming PRIM_INT__d_s_s_s[8 * 6];
273                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
274                     {
275 
276                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
277                         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]);
278                         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]);
279 
280                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
281                         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]);
282 
283                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
284                         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]);
285 
286                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
287                         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]);
288                         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]);
289 
290                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
291                         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]);
292 
293                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
294                         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]);
295                         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]);
296 
297                     }
298 
299 
300 
301                     // Forming PRIM_INT__f_s_s_s[7 * 10];
302                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
303                     {
304 
305                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
306                         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]);
307                         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]);
308 
309                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
310                         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]);
311 
312                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
313                         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]);
314 
315                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
316                         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]);
317 
318                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
319                         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]);
320 
321                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
322                         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]);
323 
324                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
325                         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]);
326                         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]);
327 
328                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
329                         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]);
330 
331                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
332                         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]);
333 
334                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
335                         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]);
336                         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]);
337 
338                     }
339 
340 
341                     VRR_K_f_s_p_s(
342                             PRIM_INT__f_s_p_s,
343                             PRIM_INT__f_s_s_s,
344                             PRIM_INT__d_s_s_s,
345                             Q_PA,
346                             aoq_PQ,
347                             one_over_2pq,
348                             3);
349 
350 
351 
352                     // Forming PRIM_INT__d_s_p_s[3 * 18];
353                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
354                     {
355 
356                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
357                         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]);
358                         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]);
359 
360                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
361                         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]);
362 
363                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
364                         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]);
365 
366                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
367                         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]);
368                         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]);
369 
370                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
371                         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]);
372                         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]);
373 
374                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
375                         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]);
376 
377                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
378                         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]);
379                         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]);
380 
381                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
382                         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]);
383 
384                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
385                         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]);
386                         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]);
387 
388                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
389                         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]);
390 
391                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
392                         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]);
393                         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]);
394 
395                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
396                         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]);
397 
398                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
399                         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]);
400 
401                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
402                         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]);
403                         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]);
404 
405                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
406                         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]);
407                         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]);
408 
409                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
410                         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]);
411 
412                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
413                         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]);
414 
415                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
416                         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]);
417                         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]);
418 
419                     }
420 
421 
422                     VRR_K_f_s_d_s(
423                             PRIM_INT__f_s_d_s,
424                             PRIM_INT__f_s_p_s,
425                             PRIM_INT__f_s_s_s,
426                             PRIM_INT__d_s_p_s,
427                             Q_PA,
428                             a_over_q,
429                             aoq_PQ,
430                             one_over_2pq,
431                             one_over_2q,
432                             2);
433 
434 
435 
436                     // Forming PRIM_INT__p_s_p_s[3 * 9];
437                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
438                     {
439 
440                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
441                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
442                         PRIM_INT__p_s_p_s[n * 9 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_p_s[n * 9 + 0]);
443 
444                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
445                         PRIM_INT__p_s_p_s[n * 9 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 1]);
446 
447                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
448                         PRIM_INT__p_s_p_s[n * 9 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_p_s[n * 9 + 2]);
449 
450                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 1]);
451                         PRIM_INT__p_s_p_s[n * 9 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 3]);
452 
453                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
454                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 4]);
455                         PRIM_INT__p_s_p_s[n * 9 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_p_s[n * 9 + 4]);
456 
457                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
458                         PRIM_INT__p_s_p_s[n * 9 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_p_s[n * 9 + 5]);
459 
460                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__p_s_s_s[n * 3 + 2]);
461                         PRIM_INT__p_s_p_s[n * 9 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 6]);
462 
463                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__p_s_s_s[n * 3 + 2]);
464                         PRIM_INT__p_s_p_s[n * 9 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 7]);
465 
466                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
467                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_p_s[n * 9 + 8]);
468                         PRIM_INT__p_s_p_s[n * 9 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_p_s[n * 9 + 8]);
469 
470                     }
471 
472 
473                     VRR_K_d_s_d_s(
474                             PRIM_INT__d_s_d_s,
475                             PRIM_INT__d_s_p_s,
476                             PRIM_INT__d_s_s_s,
477                             PRIM_INT__p_s_p_s,
478                             Q_PA,
479                             a_over_q,
480                             aoq_PQ,
481                             one_over_2pq,
482                             one_over_2q,
483                             2);
484 
485 
486                     ostei_general_vrr_K(3, 0, 3, 0, 1,
487                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
488                             PRIM_INT__f_s_d_s, PRIM_INT__f_s_p_s, NULL, PRIM_INT__d_s_d_s, NULL, PRIM_INT__f_s_f_s);
489 
490 
491                     VRR_I_g_s_s_s(
492                             PRIM_INT__g_s_s_s,
493                             PRIM_INT__f_s_s_s,
494                             PRIM_INT__d_s_s_s,
495                             P_PA,
496                             a_over_p,
497                             aop_PQ,
498                             one_over_2p,
499                             6);
500 
501 
502                     VRR_K_g_s_p_s(
503                             PRIM_INT__g_s_p_s,
504                             PRIM_INT__g_s_s_s,
505                             PRIM_INT__f_s_s_s,
506                             Q_PA,
507                             aoq_PQ,
508                             one_over_2pq,
509                             3);
510 
511 
512                     ostei_general_vrr_K(4, 0, 2, 0, 2,
513                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
514                             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);
515 
516 
517                     ostei_general_vrr_K(4, 0, 3, 0, 1,
518                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
519                             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);
520 
521 
522                     VRR_I_h_s_s_s(
523                             PRIM_INT__h_s_s_s,
524                             PRIM_INT__g_s_s_s,
525                             PRIM_INT__f_s_s_s,
526                             P_PA,
527                             a_over_p,
528                             aop_PQ,
529                             one_over_2p,
530                             5);
531 
532 
533                     ostei_general_vrr_K(5, 0, 1, 0, 3,
534                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
535                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
536 
537 
538                     ostei_general_vrr_K(5, 0, 2, 0, 2,
539                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
540                             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);
541 
542 
543                     ostei_general_vrr_K(5, 0, 3, 0, 1,
544                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
545                             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);
546 
547 
548                     ostei_general_vrr1_I(6, 4,
549                             one_over_2p, a_over_p, aop_PQ, P_PA,
550                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
551 
552 
553                     ostei_general_vrr_K(6, 0, 1, 0, 3,
554                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
555                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
556 
557 
558                     ostei_general_vrr_K(6, 0, 2, 0, 2,
559                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
560                             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);
561 
562 
563                     ostei_general_vrr_K(6, 0, 3, 0, 1,
564                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
565                             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);
566 
567 
568 
569 
570                     ////////////////////////////////////
571                     // Accumulate contracted integrals
572                     ////////////////////////////////////
573                     if(lastoffset == 0)
574                     {
575                         contract_all(100, PRIM_INT__f_s_f_s, PRIM_PTR_INT__f_s_f_s);
576                         contract_all(150, PRIM_INT__g_s_f_s, PRIM_PTR_INT__g_s_f_s);
577                         contract_all(210, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
578                         contract_all(280, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
579                     }
580                     else
581                     {
582                         contract(100, shelloffsets, PRIM_INT__f_s_f_s, PRIM_PTR_INT__f_s_f_s);
583                         contract(150, shelloffsets, PRIM_INT__g_s_f_s, PRIM_PTR_INT__g_s_f_s);
584                         contract(210, shelloffsets, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
585                         contract(280, shelloffsets, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
586                         PRIM_PTR_INT__f_s_f_s += lastoffset*100;
587                         PRIM_PTR_INT__g_s_f_s += lastoffset*150;
588                         PRIM_PTR_INT__h_s_f_s += lastoffset*210;
589                         PRIM_PTR_INT__i_s_f_s += lastoffset*280;
590                     }
591 
592                 }  // close loop over j
593             }  // close loop over i
594 
595             //Advance to the next batch
596             jstart = SIMINT_SIMD_ROUND(jend);
597 
598             //////////////////////////////////////////////
599             // Contracted integrals: Horizontal recurrance
600             //////////////////////////////////////////////
601 
602 
603             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
604 
605 
606             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
607             {
608 
609                 // set up HRR pointers
610                 double const * restrict HRR_INT__f_s_f_s = INT__f_s_f_s + abcd * 100;
611                 double const * restrict HRR_INT__g_s_f_s = INT__g_s_f_s + abcd * 150;
612                 double const * restrict HRR_INT__h_s_f_s = INT__h_s_f_s + abcd * 210;
613                 double const * restrict HRR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
614                 double * restrict HRR_INT__f_f_f_s = INT__f_f_f_s + real_abcd * 1000;
615 
616                 // form INT__f_p_f_s
617                 HRR_J_f_p(
618                     HRR_INT__f_p_f_s,
619                     HRR_INT__f_s_f_s,
620                     HRR_INT__g_s_f_s,
621                     hAB, 10);
622 
623                 // form INT__g_p_f_s
624                 HRR_J_g_p(
625                     HRR_INT__g_p_f_s,
626                     HRR_INT__g_s_f_s,
627                     HRR_INT__h_s_f_s,
628                     hAB, 10);
629 
630                 // form INT__h_p_f_s
631                 ostei_general_hrr_J(5, 1, 3, 0, hAB, HRR_INT__i_s_f_s, HRR_INT__h_s_f_s, HRR_INT__h_p_f_s);
632 
633                 // form INT__f_d_f_s
634                 HRR_J_f_d(
635                     HRR_INT__f_d_f_s,
636                     HRR_INT__f_p_f_s,
637                     HRR_INT__g_p_f_s,
638                     hAB, 10);
639 
640                 // form INT__g_d_f_s
641                 ostei_general_hrr_J(4, 2, 3, 0, hAB, HRR_INT__h_p_f_s, HRR_INT__g_p_f_s, HRR_INT__g_d_f_s);
642 
643                 // form INT__f_f_f_s
644                 ostei_general_hrr_J(3, 3, 3, 0, hAB, HRR_INT__g_d_f_s, HRR_INT__f_d_f_s, HRR_INT__f_f_f_s);
645 
646 
647             }  // close HRR loop
648 
649 
650         }   // close loop cdbatch
651 
652         istart = iend;
653     }  // close loop over ab
654 
655     return P.nshell12_clip * Q.nshell12_clip;
656 }
657 
658