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_p_f_s_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__p_f_s_p)8 int ostei_p_f_s_p(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__p_f_s_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__p_f_s_p);
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__s_f_s_p = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__s_g_s_p = work + (SIMINT_NSHELL_SIMD * 30);
30     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*75);
31     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
32     SIMINT_DBLTYPE * const restrict PRIM_INT__s_p_s_s = primwork + 6;
33     SIMINT_DBLTYPE * const restrict PRIM_INT__s_d_s_s = primwork + 21;
34     SIMINT_DBLTYPE * const restrict PRIM_INT__s_f_s_s = primwork + 45;
35     SIMINT_DBLTYPE * const restrict PRIM_INT__s_f_s_p = primwork + 75;
36     SIMINT_DBLTYPE * const restrict PRIM_INT__s_g_s_s = primwork + 105;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__s_g_s_p = primwork + 135;
38     double * const hrrwork = (double *)(primwork + 180);
39 
40 
41     // Create constants
42     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
43     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
44     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
45     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
46     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
47 
48 
49     ////////////////////////////////////////
50     // Loop over shells and primitives
51     ////////////////////////////////////////
52 
53     real_abcd = 0;
54     istart = 0;
55     for(ab = 0; ab < P.nshell12_clip; ++ab)
56     {
57         const int iend = istart + P.nprim12[ab];
58 
59         cd = 0;
60         jstart = 0;
61 
62         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
63         {
64             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
65             int jend = jstart;
66             for(i = 0; i < nshellbatch; i++)
67                 jend += Q.nprim12[cd+i];
68 
69             // Clear the beginning of the workspace (where we are accumulating integrals)
70             memset(work, 0, SIMINT_NSHELL_SIMD * 75 * sizeof(double));
71             abcd = 0;
72 
73 
74             for(i = istart; i < iend; ++i)
75             {
76                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
77 
78                 if(check_screen)
79                 {
80                     // Skip this whole thing if always insignificant
81                     if((P.screen[i] * Q.screen_max) < screen_tol)
82                         continue;
83                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
84                 }
85 
86                 icd = 0;
87                 iprimcd = 0;
88                 nprim_icd = Q.nprim12[cd];
89                 double * restrict PRIM_PTR_INT__s_f_s_p = INT__s_f_s_p + abcd * 30;
90                 double * restrict PRIM_PTR_INT__s_g_s_p = INT__s_g_s_p + abcd * 45;
91 
92 
93 
94                 // Load these one per loop over i
95                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
96                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
97                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
98 
99                 const SIMINT_DBLTYPE P_PB[3] = { SIMINT_DBLSET1(P.PB_x[i]), SIMINT_DBLSET1(P.PB_y[i]), SIMINT_DBLSET1(P.PB_z[i]) };
100 
101                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
102                 {
103                     // calculate the shell offsets
104                     // these are the offset from the shell pointed to by cd
105                     // for each element
106                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
107                     int lastoffset = 0;
108                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
109 
110                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
111                     {
112                         // Handle if the first element of the vector is a new shell
113                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
114                         {
115                             nprim_icd += Q.nprim12[cd + (++icd)];
116                             PRIM_PTR_INT__s_f_s_p += 30;
117                             PRIM_PTR_INT__s_g_s_p += 45;
118                         }
119                         iprimcd++;
120                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
121                         {
122                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
123                             {
124                                 shelloffsets[n] = shelloffsets[n-1] + 1;
125                                 lastoffset++;
126                                 nprim_icd += Q.nprim12[cd + (++icd)];
127                             }
128                             else
129                                 shelloffsets[n] = shelloffsets[n-1];
130                             iprimcd++;
131                         }
132                     }
133                     else
134                         iprimcd += SIMINT_SIMD_LEN;
135 
136                     // Do we have to compute this vector (or has it been screened out)?
137                     // (not_screened != 0 means we have to do this vector)
138                     if(check_screen)
139                     {
140                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
141                         if(vmax < screen_tol)
142                         {
143                             PRIM_PTR_INT__s_f_s_p += lastoffset*30;
144                             PRIM_PTR_INT__s_g_s_p += lastoffset*45;
145                             continue;
146                         }
147                     }
148 
149                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
150                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
151                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
152                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
153 
154 
155                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
156                     SIMINT_DBLTYPE PQ[3];
157                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
158                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
159                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
160                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
161                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
162                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
163 
164                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
165                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
166                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
167                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
168                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
169                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
170                     const SIMINT_DBLTYPE Q_PB[3] = { SIMINT_DBLLOAD(Q.PB_x, j), SIMINT_DBLLOAD(Q.PB_y, j), SIMINT_DBLLOAD(Q.PB_z, j) };
171 
172                     // NOTE: Minus sign!
173                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
174                     SIMINT_DBLTYPE aop_PQ[3];
175                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
176                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
177                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
178 
179                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
180                     SIMINT_DBLTYPE aoq_PQ[3];
181                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
182                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
183                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
184                     // Put a minus sign here so we don't have to in RR routines
185                     a_over_q = SIMINT_NEG(a_over_q);
186 
187 
188                     //////////////////////////////////////////////
189                     // Fjt function section
190                     // Maximum v value: 5
191                     //////////////////////////////////////////////
192                     // The parameter to the Fjt function
193                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
194 
195 
196                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
197 
198 
199                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 5);
200                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
201                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
202                     for(n = 0; n <= 5; n++)
203                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
204 
205                     //////////////////////////////////////////////
206                     // Primitive integrals: Vertical recurrance
207                     //////////////////////////////////////////////
208 
209                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
210                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
211                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
212                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
213                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
214                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
215                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
216 
217 
218 
219                     // Forming PRIM_INT__s_p_s_s[5 * 3];
220                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
221                     {
222 
223                         PRIM_INT__s_p_s_s[n * 3 + 0] = SIMINT_MUL(P_PB[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
224                         PRIM_INT__s_p_s_s[n * 3 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_s_s[n * 3 + 0]);
225 
226                         PRIM_INT__s_p_s_s[n * 3 + 1] = SIMINT_MUL(P_PB[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
227                         PRIM_INT__s_p_s_s[n * 3 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_s_s[n * 3 + 1]);
228 
229                         PRIM_INT__s_p_s_s[n * 3 + 2] = SIMINT_MUL(P_PB[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
230                         PRIM_INT__s_p_s_s[n * 3 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_p_s_s[n * 3 + 2]);
231 
232                     }
233 
234 
235 
236                     // Forming PRIM_INT__s_d_s_s[4 * 6];
237                     for(n = 0; n < 4; ++n)  // loop over orders of auxiliary function
238                     {
239 
240                         PRIM_INT__s_d_s_s[n * 6 + 0] = SIMINT_MUL(P_PB[0], PRIM_INT__s_p_s_s[n * 3 + 0]);
241                         PRIM_INT__s_d_s_s[n * 6 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_s_s[n * 6 + 0]);
242                         PRIM_INT__s_d_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__s_d_s_s[n * 6 + 0]);
243 
244                         PRIM_INT__s_d_s_s[n * 6 + 1] = SIMINT_MUL(P_PB[1], PRIM_INT__s_p_s_s[n * 3 + 0]);
245                         PRIM_INT__s_d_s_s[n * 6 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_s_s[n * 6 + 1]);
246 
247                         PRIM_INT__s_d_s_s[n * 6 + 2] = SIMINT_MUL(P_PB[2], PRIM_INT__s_p_s_s[n * 3 + 0]);
248                         PRIM_INT__s_d_s_s[n * 6 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_d_s_s[n * 6 + 2]);
249 
250                         PRIM_INT__s_d_s_s[n * 6 + 3] = SIMINT_MUL(P_PB[1], PRIM_INT__s_p_s_s[n * 3 + 1]);
251                         PRIM_INT__s_d_s_s[n * 6 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_d_s_s[n * 6 + 3]);
252                         PRIM_INT__s_d_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__s_d_s_s[n * 6 + 3]);
253 
254                         PRIM_INT__s_d_s_s[n * 6 + 4] = SIMINT_MUL(P_PB[2], PRIM_INT__s_p_s_s[n * 3 + 1]);
255                         PRIM_INT__s_d_s_s[n * 6 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_d_s_s[n * 6 + 4]);
256 
257                         PRIM_INT__s_d_s_s[n * 6 + 5] = SIMINT_MUL(P_PB[2], PRIM_INT__s_p_s_s[n * 3 + 2]);
258                         PRIM_INT__s_d_s_s[n * 6 + 5] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_d_s_s[n * 6 + 5]);
259                         PRIM_INT__s_d_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__s_d_s_s[n * 6 + 5]);
260 
261                     }
262 
263 
264 
265                     // Forming PRIM_INT__s_f_s_s[3 * 10];
266                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
267                     {
268 
269                         PRIM_INT__s_f_s_s[n * 10 + 0] = SIMINT_MUL(P_PB[0], PRIM_INT__s_d_s_s[n * 6 + 0]);
270                         PRIM_INT__s_f_s_s[n * 10 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_f_s_s[n * 10 + 0]);
271                         PRIM_INT__s_f_s_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_p_s_s[(n+1) * 3 + 0], PRIM_INT__s_p_s_s[n * 3 + 0]), PRIM_INT__s_f_s_s[n * 10 + 0]);
272 
273                         PRIM_INT__s_f_s_s[n * 10 + 1] = SIMINT_MUL(P_PB[1], PRIM_INT__s_d_s_s[n * 6 + 0]);
274                         PRIM_INT__s_f_s_s[n * 10 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_f_s_s[n * 10 + 1]);
275 
276                         PRIM_INT__s_f_s_s[n * 10 + 2] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 0]);
277                         PRIM_INT__s_f_s_s[n * 10 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 0], PRIM_INT__s_f_s_s[n * 10 + 2]);
278 
279                         PRIM_INT__s_f_s_s[n * 10 + 3] = SIMINT_MUL(P_PB[0], PRIM_INT__s_d_s_s[n * 6 + 3]);
280                         PRIM_INT__s_f_s_s[n * 10 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_f_s_s[n * 10 + 3]);
281 
282                         PRIM_INT__s_f_s_s[n * 10 + 4] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 1]);
283                         PRIM_INT__s_f_s_s[n * 10 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 1], PRIM_INT__s_f_s_s[n * 10 + 4]);
284 
285                         PRIM_INT__s_f_s_s[n * 10 + 5] = SIMINT_MUL(P_PB[0], PRIM_INT__s_d_s_s[n * 6 + 5]);
286                         PRIM_INT__s_f_s_s[n * 10 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_f_s_s[n * 10 + 5]);
287 
288                         PRIM_INT__s_f_s_s[n * 10 + 6] = SIMINT_MUL(P_PB[1], PRIM_INT__s_d_s_s[n * 6 + 3]);
289                         PRIM_INT__s_f_s_s[n * 10 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_f_s_s[n * 10 + 6]);
290                         PRIM_INT__s_f_s_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_p_s_s[(n+1) * 3 + 1], PRIM_INT__s_p_s_s[n * 3 + 1]), PRIM_INT__s_f_s_s[n * 10 + 6]);
291 
292                         PRIM_INT__s_f_s_s[n * 10 + 7] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 3]);
293                         PRIM_INT__s_f_s_s[n * 10 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 3], PRIM_INT__s_f_s_s[n * 10 + 7]);
294 
295                         PRIM_INT__s_f_s_s[n * 10 + 8] = SIMINT_MUL(P_PB[1], PRIM_INT__s_d_s_s[n * 6 + 5]);
296                         PRIM_INT__s_f_s_s[n * 10 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_f_s_s[n * 10 + 8]);
297 
298                         PRIM_INT__s_f_s_s[n * 10 + 9] = SIMINT_MUL(P_PB[2], PRIM_INT__s_d_s_s[n * 6 + 5]);
299                         PRIM_INT__s_f_s_s[n * 10 + 9] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_d_s_s[(n+1) * 6 + 5], PRIM_INT__s_f_s_s[n * 10 + 9]);
300                         PRIM_INT__s_f_s_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_p_s_s[(n+1) * 3 + 2], PRIM_INT__s_p_s_s[n * 3 + 2]), PRIM_INT__s_f_s_s[n * 10 + 9]);
301 
302                     }
303 
304 
305                     VRR_L_s_f_s_p(
306                             PRIM_INT__s_f_s_p,
307                             PRIM_INT__s_f_s_s,
308                             PRIM_INT__s_d_s_s,
309                             Q_PB,
310                             aoq_PQ,
311                             one_over_2pq,
312                             1);
313 
314 
315                     VRR_J_s_g_s_s(
316                             PRIM_INT__s_g_s_s,
317                             PRIM_INT__s_f_s_s,
318                             PRIM_INT__s_d_s_s,
319                             P_PB,
320                             a_over_p,
321                             aop_PQ,
322                             one_over_2p,
323                             2);
324 
325 
326                     VRR_L_s_g_s_p(
327                             PRIM_INT__s_g_s_p,
328                             PRIM_INT__s_g_s_s,
329                             PRIM_INT__s_f_s_s,
330                             Q_PB,
331                             aoq_PQ,
332                             one_over_2pq,
333                             1);
334 
335 
336 
337 
338                     ////////////////////////////////////
339                     // Accumulate contracted integrals
340                     ////////////////////////////////////
341                     if(lastoffset == 0)
342                     {
343                         contract_all(30, PRIM_INT__s_f_s_p, PRIM_PTR_INT__s_f_s_p);
344                         contract_all(45, PRIM_INT__s_g_s_p, PRIM_PTR_INT__s_g_s_p);
345                     }
346                     else
347                     {
348                         contract(30, shelloffsets, PRIM_INT__s_f_s_p, PRIM_PTR_INT__s_f_s_p);
349                         contract(45, shelloffsets, PRIM_INT__s_g_s_p, PRIM_PTR_INT__s_g_s_p);
350                         PRIM_PTR_INT__s_f_s_p += lastoffset*30;
351                         PRIM_PTR_INT__s_g_s_p += lastoffset*45;
352                     }
353 
354                 }  // close loop over j
355             }  // close loop over i
356 
357             //Advance to the next batch
358             jstart = SIMINT_SIMD_ROUND(jend);
359 
360             //////////////////////////////////////////////
361             // Contracted integrals: Horizontal recurrance
362             //////////////////////////////////////////////
363 
364 
365             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
366 
367 
368             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
369             {
370 
371                 // set up HRR pointers
372                 double const * restrict HRR_INT__s_f_s_p = INT__s_f_s_p + abcd * 30;
373                 double const * restrict HRR_INT__s_g_s_p = INT__s_g_s_p + abcd * 45;
374                 double * restrict HRR_INT__p_f_s_p = INT__p_f_s_p + real_abcd * 90;
375 
376                 // form INT__p_f_s_p
377                 HRR_I_p_f(
378                     HRR_INT__p_f_s_p,
379                     HRR_INT__s_f_s_p,
380                     HRR_INT__s_g_s_p,
381                     hAB, 3);
382 
383 
384             }  // close HRR loop
385 
386 
387         }   // close loop cdbatch
388 
389         istart = iend;
390     }  // close loop over ab
391 
392     return P.nshell12_clip * Q.nshell12_clip;
393 }
394 
395